### The Programmatic Coup: Ideology, the Military, and Political Violence
### Version of 30 June 2025
### Kevin Koehler, kevin.koehler@santannapisa.it

rm(list=ls())

# Packages ----------------------------------------------------------------
library(foreign)
library(tidyverse)
library(readstata13)
library(readxl)
library(countrycode)
library(miceadds)
library(stargazer)
library(jtools)
library(sandwich)
library(lmtest)
library(tikzDevice)
library(rio)
library(dttr2)
library(car)

# CAM data ----------------------------------------------------------------
cam <- read.dta13("https://militarycoups.org/CAM_Ideology_v1.dta")
cam <- cam %>%
  mutate(cowcode=ccode) %>%
  select(-ccode) %>%
  relocate(cowcode) %>%
  mutate(date1 = as.Date(date1, format = "%d%b%Y")) %>%
  mutate(date2 = as.Date(date2, format = "%d%b%Y")) %>%
  mutate(date3 = as.Date(date3, format = "%d%b%Y")) %>%
  mutate(date4 = as.Date(date4, format = "%d%b%Y")) 

# Descriptive figures ------------------------------------------------------

## Figure 2 ----
region <- cam %>%
  group_by(region) %>%
  mutate(left=sum(left1, na.rm=T)+sum(left2, na.rm=T)+sum(left3, na.rm=T)+sum(left4, na.rm=T),
         right=sum(right1, na.rm=T)+sum(right2, na.rm=T)+sum(right3, na.rm=T)+sum(right4, na.rm=T),
         coup=sum(successful1, na.rm=T)+sum(successful2, na.rm=T)+sum(successful3, na.rm=T)+sum(successful4, na.rm=T),
         non_ideo=coup-left-right,
         region=factor(region,
                       levels = c(1,2,3,4,5),
                       labels = c("Europe","MENA","Sub-Saharan Africa","Asia","Americas"))) %>%
  select(region, left, right, non_ideo) %>%
  distinct(region, .keep_all = T) %>%
  pivot_longer(cols = c("left","right","non_ideo"),
               names_to = "type") 

region_fig <- ggplot(region, aes(x=reorder(region,value),
                                 y=value, 
                                 fill=type)) +
  geom_bar(stat="identity",
           position = "dodge") +
  theme_classic() +
  xlab("") +
  ylab("") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text = element_text(family="Times New Roman")) +
  scale_fill_brewer(palette = "Greys",
                    labels = c("Left-wing", 
                                 "Non-ideological", 
                                 "Right-wing")) +
  coord_flip()
region_fig
ggsave("region.jpg", dpi=300)

## Figure 3 ----
decade <- cam %>%
  mutate(decade=floor(year/10)*10) %>%
  group_by(decade) %>%
  mutate(left=sum(left1, na.rm=T)+sum(left2, na.rm=T)+sum(left3, na.rm=T)+sum(left4, na.rm=T),
         right=sum(right1, na.rm=T)+sum(right2, na.rm=T)+sum(right3, na.rm=T)+sum(right4, na.rm=T),
         coup=sum(successful1, na.rm=T)+sum(successful2, na.rm=T)+sum(successful3, na.rm=T)+sum(successful4, na.rm=T),
         non_ideo=coup-left-right) %>%
  select(decade, left, right, non_ideo) %>%
  distinct(decade, .keep_all = T) %>%
  pivot_longer(cols = c("left","right","non_ideo"),
               names_to = "type") %>%
  mutate(decade=factor(decade,
                       levels = c(1950,1960,1970,1980,1990,2000,2010,2020),
                       labels = c("1950s","1960s","1970s","1980s","1990s","2000s","2010s","2020s")))

decade_fig <- ggplot(decade, aes(x=decade,
                                 y=value, 
                                 fill=type)) +
  geom_bar(stat="identity") +
  scale_x_discrete(limits=rev) +
  theme_classic() +
  xlab("") +
  ylab("") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text = element_text(family="Times New Roman")) +
  scale_fill_brewer(palette = "Greys",
                    breaks = c("left", 
                               "non_ideo", 
                               "right"),
                    labels = c("Left-wing", 
                               "Non-ideological", 
                               "Right-wing")) +
  coord_flip()
decade_fig
ggsave("decade.jpg", dpi=300)

# Overlap between CAM ideology and COLPUS types ----

list <- read.csv("https://www.militarycoups.org/cam_list_current.csv")
list$left <- NA_real_
list$right <- NA_real_
for(c in 1:length(row_number(list))){
  for(d in 1:length(row_number(cam))){
    if(!is.na(cam$date1[d])) {
      if(list$cowcode[c]==cam$cowcode[d] & list$date[c]==cam$date1[d]) {
        list$left[c] <- cam$left1[d]
        list$right[c] <- cam$right1[d]
      }
      if(!is.na(cam$date2[d]) & list$cowcode[c]==cam$cowcode[d] & list$date[c]==cam$date2[d]) {
        list$left[c] <- cam$left2[d]
        list$right[c] <- cam$right2[d]
      }
      if(!is.na(cam$date3[d]) & list$cowcode[c]==cam$cowcode[d] & list$date[c]==cam$date3[d]) {
        list$left[c] <- cam$left3[d]
        list$right[c] <- cam$right3[d]
      }
      if(!is.na(cam$date4[d]) & list$cowcode[c]==cam$cowcode[d] & list$date[c]==cam$date4[d]) {
        list$left[c] <- cam$left4[d]
        list$right[c] <- cam$right4[d]
      }
    }
  }
}

## Check if merging worked OK ----

cam$right_sum <- rowSums(cam[, c("right1", "right2", "right3", "right4")], na.rm = T)
sum(cam$right_sum)
table(list$right)

## same number of right-wing coups

cam$left_sum <- rowSums(cam[, c("left1", "left2", "left3", "left4")], na.rm = T)
sum(cam$left_sum)
table(list$left)

## same number of left-wing coups

list <- list %>%
  mutate(date=as.Date(date)) %>%
  select(-coup) %>%
  mutate(cam=1)

## Colpus ----
### data from https://doi.org/10.1093/isq/sqab058
colpus <- import("/Users/kevin/Dropbox/CAM Dataset/Other Datasets/COLPUS/Data and Code/Master_Coup List_all regimes_basic.dta")
colpus <- colpus %>%
  rename(cowcode=ccode) %>%
  filter(milcoup==1) %>%
  select(-coup) %>%
  mutate(colpus=1)

comp <- full_join(list, colpus, by=c("cowcode","date","year"))

comp <- comp %>%
  mutate(colpus=ifelse(is.na(colpus),0,colpus),
         cam=ifelse(is.na(cam),0,cam),
         p_coup=case_when(
           left==1 | right==1~1,
           left==0 | right==0~0,
           TRUE~NA_real_
         ))

## This is the information in footnote 1:
table(comp$type2,comp$p_coup)

# IVs based on CAM ----

## Ideology of last successful coup in year ----
ideo <- cam %>%
  mutate(any_coup=case_when(
    rowSums(cbind(successful1,successful2,successful3,successful4),na.rm = T)>0~1,
    TRUE~0
  )) %>%
  mutate(last_coup=rowSums(cbind(successful1,successful2,successful3,successful4),na.rm = T)) %>%
  mutate(ideo_last_coup=case_when(
    last_coup==1 & right1==1~"right",
    last_coup==2 & right2==1~"right",
    last_coup==3 & right3==1~"right",
    last_coup==4 & right4==1~"right",
    last_coup==1 & left1==1~"left",
    last_coup==2 & left2==1~"left",
    last_coup==3 & left3==1~"left",
    last_coup==4 & left4==1~"left",
    TRUE~NA_character_
  )) %>%
  mutate(date=case_when(
    last_coup==1~date1,
    last_coup==2~date2,
    last_coup==3~date3,
    last_coup==4~date4,
    TRUE~NA_Date_
  )) %>%
  mutate(agency_last_coup=case_when(
    last_coup==1 & junior1==1~1,
    last_coup==2 & junior2==1~1,
    last_coup==3 & junior3==1~1,
    last_coup==4 & junior4==1~1,
    last_coup==1 & junior1==0~0,
    last_coup==2 & junior2==0~0,
    last_coup==3 & junior3==0~0,
    last_coup==4 & junior4==0~0,
    TRUE~NA_real_
  )) %>%
  mutate(casualties=case_when(
    last_coup==1~casult1,
    last_coup==2~casult2,
    last_coup==3~casult3,
    last_coup==4~casult4,
    TRUE~NA_integer_
  )) %>%
  select(cowcode, year, any_coup, last_coup, agency_last_coup, ideo_last_coup, date, casualties) %>%
  mutate(i_coup=ifelse(!is.na(ideo_last_coup),1,0)) %>%
  filter(any_coup==1) 

## Three-year windows after successful coups ----
list <- read.csv("https://www.militarycoups.org/cam_list_current.csv") %>%
  filter(successful==1) %>% ##Only successful coups
  group_by(cowcode, year) %>% 
  filter(no==max(no)) ## Only last successful coup if several coups in same year

episodes <- data.frame()
for(i in 1:length(list$cowcode)) {
  id <- rep(paste(list$country[i],list$year[i],sep="_"),4)
  cowcode <- rep(list$cowcode[i],4)
  country <- rep(list$country[i],4)
  combat <- c(list$combat[i],NA,NA,NA)
  coup <- c(1,0,0,0)
  coup_year <- rep(list$year[i],4)
  year <- c(list$year[i],
            list$year[i]+1,
            list$year[i]+2,
            list$year[i]+3)
  episode <- data.frame(id, cowcode, country, year, coup, combat, coup_year)
  episodes <- rbind(episodes, episode)
}

### flag overlapping episodes ----
all_episodes <- paste(episodes$country,episodes$year,sep="_")
for(i in 1:length(episodes$id)) {
  episodes$overlap[i] <- 0
  for(n in 1:3) {
    if(table(all_episodes)[names(table(all_episodes))==paste(episodes$country[i],episodes$coup_year[i]+n,sep="_")]>1) {
      episodes$overlap[i] <- 1 
    }
  }
}
episodes <- episodes %>%
  filter(overlap==0) %>%
  group_by(id) %>%
  mutate(y_post_coup=sequence(n())-1)

episodes <- left_join(episodes, ideo, by=c("cowcode","year"))

### Constructing IVs ----

episodes <- episodes %>%
  mutate(i_coup=ifelse(coup==1 & is.na(i_coup),0,i_coup)) %>%
  group_by(id) %>%
  mutate(p_coup=case_when(
    sum(i_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(any_coup=case_when(
    sum(any_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(coup_agency=case_when(
    sum(combat==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(casualties=as.numeric(casualties)) %>%
  mutate(casualties=case_when(
    y_post_coup>0~sum(casualties,na.rm = T),
    TRUE~0
  )) %>%
  select(cowcode, year, p_coup, any_coup, id, coup, combat, coup_agency, casualties, y_post_coup)

## Two-year windows after successful coups ----

list <- read.csv("https://www.militarycoups.org/cam_list_current.csv") %>%
  filter(successful==1) %>% ##Only successful coups
  group_by(cowcode, year) %>% 
  filter(no==max(no)) ## Only last successful coup if several coups in same year

episodes_2 <- data.frame()
for(i in 1:length(list$cowcode)) {
  id <- rep(paste(list$country[i],list$year[i],sep="_"),3)
  cowcode <- rep(list$cowcode[i],3)
  country <- rep(list$country[i],3)
  combat <- c(list$combat[i],NA,NA)
  coup <- c(1,0,0)
  coup_year <- rep(list$year[i],3)
  year <- c(list$year[i],
            list$year[i]+1,
            list$year[i]+2)
  episode <- data.frame(id, cowcode, country, year, coup, combat, coup_year)
  episodes_2 <- rbind(episodes_2, episode)
}

### flag overlapping episodes ----
all_episodes_2 <- paste(episodes_2$country,episodes_2$year,sep="_")
for(i in 1:length(episodes_2$id)) {
  episodes_2$overlap[i] <- 0
  for(n in 1:2) {
    if(table(all_episodes_2)[names(table(all_episodes_2))==paste(episodes_2$country[i],episodes_2$coup_year[i]+n,sep="_")]>1) {
      episodes_2$overlap[i] <- 1 
    }
  }
}
episodes_2 <- episodes_2 %>%
  filter(overlap==0) %>%
  group_by(id) %>%
  mutate(y_post_coup=sequence(n())-1)

episodes_2 <- left_join(episodes_2, ideo, by=c("cowcode","year"))

### Constructing IVs ----

episodes_2 <- episodes_2 %>%
  mutate(i_coup=ifelse(coup==1 & is.na(i_coup),0,i_coup)) %>%
  group_by(id) %>%
  mutate(p_coup=case_when(
    sum(i_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(any_coup=case_when(
    sum(any_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(coup_agency=case_when(
    sum(combat==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(casualties=as.numeric(casualties)) %>%
  mutate(casualties=case_when(
    y_post_coup>0~sum(casualties,na.rm = T),
    TRUE~0
  )) %>%
  select(cowcode, year, p_coup, any_coup, id, coup, combat, coup_agency, casualties, y_post_coup)

## One-year windows after successful coups ----

list <- read.csv("https://www.militarycoups.org/cam_list_current.csv") %>%
  filter(successful==1) %>% ##Only successful coups
  group_by(cowcode, year) %>% 
  filter(no==max(no)) ## Only last successful coup if several coups in same year

episodes_1 <- data.frame()
for(i in 1:length(list$cowcode)) {
  id <- rep(paste(list$country[i],list$year[i],sep="_"),2)
  cowcode <- rep(list$cowcode[i],2)
  country <- rep(list$country[i],2)
  combat <- c(list$combat[i],NA)
  coup <- c(1,0)
  coup_year <- rep(list$year[i],2)
  year <- c(list$year[i],
            list$year[i]+1)
  episode <- data.frame(id, cowcode, country, year, coup, combat, coup_year)
  episodes_1 <- rbind(episodes_1, episode)
}

### flag overlapping episodes ----
all_episodes_1 <- paste(episodes_1$country,episodes_1$year,sep="_")
for(i in 1:length(episodes_1$id)) {
  episodes_1$overlap[i] <- 0
  for(n in 1:1) {
    if(table(all_episodes_1)[names(table(all_episodes_1))==paste(episodes_1$country[i],episodes_1$coup_year[i]+n,sep="_")]>1) {
      episodes_1$overlap[i] <- 1 
    }
  }
}
episodes_1 <- episodes_1 %>%
  filter(overlap==0) %>%
  group_by(id) %>%
  mutate(y_post_coup=sequence(n())-1)

episodes_1 <- left_join(episodes_1, ideo, by=c("cowcode","year"))

### Constructing IVs ----

episodes_1 <- episodes_1 %>%
  mutate(i_coup=ifelse(coup==1 & is.na(i_coup),0,i_coup)) %>%
  group_by(id) %>%
  mutate(p_coup=case_when(
    sum(i_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(any_coup=case_when(
    sum(any_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(coup_agency=case_when(
    sum(combat==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(casualties=as.numeric(casualties)) %>%
  mutate(casualties=case_when(
    y_post_coup>0~sum(casualties,na.rm = T),
    TRUE~0
  )) %>%
  select(cowcode, year, p_coup, any_coup, id, coup, combat, coup_agency, casualties, y_post_coup)

## Four-year windows after successful coups ----

list <- read.csv("https://www.militarycoups.org/cam_list_current.csv") %>%
  filter(successful==1) %>% ##Only successful coups
  group_by(cowcode, year) %>% 
  filter(no==max(no)) ## Only last successful coup if several coups in same year

episodes_4 <- data.frame()
for(i in 1:length(list$cowcode)) {
  id <- rep(paste(list$country[i],list$year[i],sep="_"),5)
  cowcode <- rep(list$cowcode[i],5)
  country <- rep(list$country[i],5)
  combat <- c(list$combat[i],NA,NA,NA,NA)
  coup <- c(1,0,0,0,0)
  coup_year <- rep(list$year[i],5)
  year <- c(list$year[i],
            list$year[i]+1,
            list$year[i]+2,
            list$year[i]+3,
            list$year[i]+4)
  episode <- data.frame(id, cowcode, country, year, coup, combat, coup_year)
  episodes_4 <- rbind(episodes_4, episode)
}

### flag overlapping episodes ----
all_episodes_4 <- paste(episodes_4$country,episodes_4$year,sep="_")
for(i in 1:length(episodes_4$id)) {
  episodes_4$overlap[i] <- 0
  for(n in 1:4) {
    if(table(all_episodes_4)[names(table(all_episodes_4))==paste(episodes_4$country[i],episodes_4$coup_year[i]+n,sep="_")]>1) {
      episodes_4$overlap[i] <- 1 
    }
  }
}
episodes_4 <- episodes_4 %>%
  filter(overlap==0) %>%
  group_by(id) %>%
  mutate(y_post_coup=sequence(n())-1)

episodes_4 <- left_join(episodes_4, ideo, by=c("cowcode","year"))

### Constructing IVs ----

episodes_4 <- episodes_4 %>%
  mutate(i_coup=ifelse(coup==1 & is.na(i_coup),0,i_coup)) %>%
  group_by(id) %>%
  mutate(p_coup=case_when(
    sum(i_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(any_coup=case_when(
    sum(any_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(coup_agency=case_when(
    sum(combat==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(casualties=as.numeric(casualties)) %>%
  mutate(casualties=case_when(
    y_post_coup>0~sum(casualties,na.rm = T),
    TRUE~0
  )) %>%
  select(cowcode, year, p_coup, any_coup, id, coup, combat, coup_agency, casualties, y_post_coup)

## Five-year windows after successful coups ----

list <- read.csv("https://www.militarycoups.org/cam_list_current.csv") %>%
  filter(successful==1) %>% ##Only successful coups
  group_by(cowcode, year) %>% 
  filter(no==max(no)) ## Only last successful coup if several coups in same year

episodes_5 <- data.frame()
for(i in 1:length(list$cowcode)) {
  id <- rep(paste(list$country[i],list$year[i],sep="_"),6)
  cowcode <- rep(list$cowcode[i],6)
  country <- rep(list$country[i],6)
  combat <- c(list$combat[i],NA,NA,NA,NA,NA)
  coup <- c(1,0,0,0,0,0)
  coup_year <- rep(list$year[i],6)
  year <- c(list$year[i],
            list$year[i]+1,
            list$year[i]+2,
            list$year[i]+3,
            list$year[i]+4,
            list$year[i]+5)
  episode <- data.frame(id, cowcode, country, year, coup, combat, coup_year)
  episodes_5 <- rbind(episodes_5, episode)
}

### flag overlapping episodes ----
all_episodes_5 <- paste(episodes_5$country,episodes_5$year,sep="_")
for(i in 1:length(episodes_5$id)) {
  episodes_5$overlap[i] <- 0
  for(n in 1:5) {
    if(table(all_episodes_5)[names(table(all_episodes_5))==paste(episodes_5$country[i],episodes_5$coup_year[i]+n,sep="_")]>1) {
      episodes_5$overlap[i] <- 1 
    }
  }
}
episodes_5 <- episodes_5 %>%
  filter(overlap==0) %>%
  group_by(id) %>%
  mutate(y_post_coup=sequence(n())-1)

episodes_5 <- left_join(episodes_5, ideo, by=c("cowcode","year"))

### Constructing IVs ----

episodes_5 <- episodes_5 %>%
  mutate(i_coup=ifelse(coup==1 & is.na(i_coup),0,i_coup)) %>%
  group_by(id) %>%
  mutate(p_coup=case_when(
    sum(i_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(any_coup=case_when(
    sum(any_coup==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(coup_agency=case_when(
    sum(combat==1,na.rm = T) & y_post_coup>0~1,
    TRUE~0
  )) %>%
  mutate(casualties=as.numeric(casualties)) %>%
  mutate(casualties=case_when(
    y_post_coup>0~sum(casualties,na.rm = T),
    TRUE~0
  )) %>%
  select(cowcode, year, p_coup, any_coup, id, coup, combat, coup_agency, casualties, y_post_coup)

# DVs ----

## VDem ----
## reversing VDem measures, renaming VDem-based DVs

## Load V-Dem-CY-Full+Others-v13 from https://www.v-dem.net/data/the-v-dem-dataset/country-year-v-dem-fullothers-v13/ into an object called data

data <- data %>%
  filter(year>1949 & year<2018) %>%
  rename(cowcode=COWcode) %>%
  mutate(torture=v2cltort*-1) %>%
  mutate(killings=v2clkill*-1) %>%
  mutate(protest=v2cagenmob)
hist(data$v2cltort)
hist(data$v2clkill)

## HR protection scores ----

## Load the Latent Human Rights Protection Scores (HumanRightsProtectionScores_v4.01.csv) from https://search.dataone.org/view/sha256:e09f5137c377d84894fc166a655ff1caf099f9e69f01039c0e27d71416097976  

hrps <- hrps %>%
  rename(year=YEAR) %>%
  rename(cowcode=COW) %>%
  mutate(repression=theta_mean*-1) %>%
  rename(repression_sd=theta_sd) %>%
  select(cowcode, year, repression, repression_sd)
data <- left_join(data,hrps, by=c("cowcode","year"))
hist(data$repression)

## Merging (main data) ----
data_main <- left_join(data, episodes, by=c("cowcode","year")) %>%
  mutate(p_coup=ifelse(is.na(p_coup),0,p_coup)) %>%
  mutate(any_coup=ifelse(is.na(any_coup),0,any_coup)) %>%
  mutate(coup_agency=ifelse(is.na(coup_agency),0,coup_agency)) %>%
  mutate(casualties=ifelse(is.na(casualties),0,casualties))

### Dropping unnecessary variables ----
data_main <- data_main %>%
  select(cowcode, 
         year, 
         p_coup, 
         any_coup, 
         y_post_coup,
         id, 
         coup, 
         combat,
         coup_agency,
         casualties,
         torture,
         killings,
         protest,
         repression, 
         e_polity2,
         v2x_polyarchy,
         v2clpolcl,
         v2regimpgroup,
         v2regsupgroupssize) %>%
  mutate(coup=ifelse(is.na(coup),0,coup)) %>%
  group_by(cowcode) %>%
  mutate(coup_id=paste(cowcode,cumsum(coup),sep="_")) %>%
  mutate(p_dem=case_when(
    coup==1 & (lag(e_polity2,n=1)>6 |
                 lag(e_polity2,n=2)>6 |
                 lag(e_polity2,n=3)>6 |
                 lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  ungroup %>%
  relocate(coup_id, .after=coup) %>%
  group_by(coup_id) %>%
  mutate(tnc=row_number()) %>%
  mutate(p_dem=case_when(
    sum(p_dem>0,na.rm = T) & y_post_coup>0~1,
    TRUE~p_dem
  )) %>%
  ungroup %>%
  relocate(tnc, .after=coup_id) %>%
  group_by(cowcode) %>%
  mutate(p_dem=case_when(
    is.na(y_post_coup) & (lag(e_polity2,n=1)>6 |
                            lag(e_polity2,n=2)>6 |
                            lag(e_polity2,n=3)>6 |
                            lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  mutate(prior_dem= 
           (lag(e_polity2,n=5)+
              lag(e_polity2,n=4)+
              lag(e_polity2,n=3)+
              lag(e_polity2,n=2))/4) %>%
  ungroup %>%
  mutate(prior_dem=ifelse(prior_dem>=6,1,0))

## Merging (two-year window) ----
data_2 <- left_join(data, episodes_2, by=c("cowcode","year")) %>%
  mutate(p_coup=ifelse(is.na(p_coup),0,p_coup)) %>%
  mutate(any_coup=ifelse(is.na(any_coup),0,any_coup)) %>%
  mutate(coup_agency=ifelse(is.na(coup_agency),0,coup_agency)) %>%
  mutate(casualties=ifelse(is.na(casualties),0,casualties))

### Dropping unnecessary variables ----
data_2 <- data_2 %>%
  select(cowcode, 
         year, 
         p_coup, 
         any_coup, 
         y_post_coup,
         id, 
         coup, 
         combat,
         coup_agency,
         casualties,
         torture,
         killings,
         protest,
         repression, 
         e_polity2,
         v2x_polyarchy,
         v2clpolcl,
         v2regimpgroup,
         v2regsupgroupssize) %>%
  mutate(coup=ifelse(is.na(coup),0,coup)) %>%
  group_by(cowcode) %>%
  mutate(coup_id=paste(cowcode,cumsum(coup),sep="_")) %>%
  mutate(p_dem=case_when(
    coup==1 & (lag(e_polity2,n=1)>6 |
                 lag(e_polity2,n=2)>6 |
                 lag(e_polity2,n=3)>6 |
                 lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  ungroup %>%
  relocate(coup_id, .after=coup) %>%
  group_by(coup_id) %>%
  mutate(tnc=row_number()) %>%
  mutate(p_dem=case_when(
    sum(p_dem>0,na.rm = T) & y_post_coup>0~1,
    TRUE~p_dem
  )) %>%
  ungroup %>%
  relocate(tnc, .after=coup_id) %>%
  group_by(cowcode) %>%
  mutate(p_dem=case_when(
    is.na(y_post_coup) & (lag(e_polity2,n=1)>6 |
                            lag(e_polity2,n=2)>6 |
                            lag(e_polity2,n=3)>6 |
                            lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  mutate(prior_dem= 
           (lag(e_polity2,n=5)+
              lag(e_polity2,n=4)+
              lag(e_polity2,n=3)+
              lag(e_polity2,n=2))/4) %>%
  ungroup %>%
  mutate(prior_dem=ifelse(prior_dem>=6,1,0))

## Merging (one-year window) ----
data_1 <- left_join(data, episodes_1, by=c("cowcode","year")) %>%
  mutate(p_coup=ifelse(is.na(p_coup),0,p_coup)) %>%
  mutate(any_coup=ifelse(is.na(any_coup),0,any_coup)) %>%
  mutate(coup_agency=ifelse(is.na(coup_agency),0,coup_agency)) %>%
  mutate(casualties=ifelse(is.na(casualties),0,casualties))

### Dropping unnecessary variables ----
data_1 <- data_1 %>%
  select(cowcode, 
         year, 
         p_coup, 
         any_coup, 
         y_post_coup,
         id, 
         coup, 
         combat,
         coup_agency,
         casualties,
         torture,
         killings,
         protest,
         repression, 
         e_polity2,
         v2x_polyarchy,
         v2clpolcl,
         v2regimpgroup,
         v2regsupgroupssize) %>%
  mutate(coup=ifelse(is.na(coup),0,coup)) %>%
  group_by(cowcode) %>%
  mutate(coup_id=paste(cowcode,cumsum(coup),sep="_")) %>%
  mutate(p_dem=case_when(
    coup==1 & (lag(e_polity2,n=1)>6 |
                 lag(e_polity2,n=2)>6 |
                 lag(e_polity2,n=3)>6 |
                 lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  ungroup %>%
  relocate(coup_id, .after=coup) %>%
  group_by(coup_id) %>%
  mutate(tnc=row_number()) %>%
  mutate(p_dem=case_when(
    sum(p_dem>0,na.rm = T) & y_post_coup>0~1,
    TRUE~p_dem
  )) %>%
  ungroup %>%
  relocate(tnc, .after=coup_id) %>%
  group_by(cowcode) %>%
  mutate(p_dem=case_when(
    is.na(y_post_coup) & (lag(e_polity2,n=1)>6 |
                            lag(e_polity2,n=2)>6 |
                            lag(e_polity2,n=3)>6 |
                            lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  mutate(prior_dem= 
           (lag(e_polity2,n=5)+
              lag(e_polity2,n=4)+
              lag(e_polity2,n=3)+
              lag(e_polity2,n=2))/4) %>%
  ungroup %>%
  mutate(prior_dem=ifelse(prior_dem>=6,1,0))

## Merging (four-year window) ----
data_4 <- left_join(data, episodes_4, by=c("cowcode","year")) %>%
  mutate(p_coup=ifelse(is.na(p_coup),0,p_coup)) %>%
  mutate(any_coup=ifelse(is.na(any_coup),0,any_coup)) %>%
  mutate(coup_agency=ifelse(is.na(coup_agency),0,coup_agency)) %>%
  mutate(casualties=ifelse(is.na(casualties),0,casualties))

### Dropping unnecessary variables ----
data_4 <- data_4 %>%
  select(cowcode, 
         year, 
         p_coup, 
         any_coup, 
         y_post_coup,
         id, 
         coup, 
         combat,
         coup_agency,
         casualties,
         torture,
         killings,
         protest,
         repression, 
         e_polity2,
         v2x_polyarchy,
         v2clpolcl,
         v2regimpgroup,
         v2regsupgroupssize) %>%
  mutate(coup=ifelse(is.na(coup),0,coup)) %>%
  group_by(cowcode) %>%
  mutate(coup_id=paste(cowcode,cumsum(coup),sep="_")) %>%
  mutate(p_dem=case_when(
    coup==1 & (lag(e_polity2,n=1)>6 |
                 lag(e_polity2,n=2)>6 |
                 lag(e_polity2,n=3)>6 |
                 lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  ungroup %>%
  relocate(coup_id, .after=coup) %>%
  group_by(coup_id) %>%
  mutate(tnc=row_number()) %>%
  mutate(p_dem=case_when(
    sum(p_dem>0,na.rm = T) & y_post_coup>0~1,
    TRUE~p_dem
  )) %>%
  ungroup %>%
  relocate(tnc, .after=coup_id) %>%
  group_by(cowcode) %>%
  mutate(p_dem=case_when(
    is.na(y_post_coup) & (lag(e_polity2,n=1)>6 |
                            lag(e_polity2,n=2)>6 |
                            lag(e_polity2,n=3)>6 |
                            lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  mutate(prior_dem= 
           (lag(e_polity2,n=5)+
              lag(e_polity2,n=4)+
              lag(e_polity2,n=3)+
              lag(e_polity2,n=2))/4) %>%
  ungroup %>%
  mutate(prior_dem=ifelse(prior_dem>=6,1,0))

## Merging (five-year window) ----
data_5 <- left_join(data, episodes_5, by=c("cowcode","year")) %>%
  mutate(p_coup=ifelse(is.na(p_coup),0,p_coup)) %>%
  mutate(any_coup=ifelse(is.na(any_coup),0,any_coup)) %>%
  mutate(coup_agency=ifelse(is.na(coup_agency),0,coup_agency)) %>%
  mutate(casualties=ifelse(is.na(casualties),0,casualties))

### Dropping unnecessary variables ----
data_5 <- data_5 %>%
  select(cowcode, 
         year, 
         p_coup, 
         any_coup, 
         y_post_coup,
         id, 
         coup, 
         combat,
         coup_agency,
         casualties,
         torture,
         killings,
         protest,
         repression, 
         e_polity2,
         v2x_polyarchy,
         v2clpolcl,
         v2regimpgroup,
         v2regsupgroupssize) %>%
  mutate(coup=ifelse(is.na(coup),0,coup)) %>%
  group_by(cowcode) %>%
  mutate(coup_id=paste(cowcode,cumsum(coup),sep="_")) %>%
  mutate(p_dem=case_when(
    coup==1 & (lag(e_polity2,n=1)>6 |
                 lag(e_polity2,n=2)>6 |
                 lag(e_polity2,n=3)>6 |
                 lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  ungroup %>%
  relocate(coup_id, .after=coup) %>%
  group_by(coup_id) %>%
  mutate(tnc=row_number()) %>%
  mutate(p_dem=case_when(
    sum(p_dem>0,na.rm = T) & y_post_coup>0~1,
    TRUE~p_dem
  )) %>%
  ungroup %>%
  relocate(tnc, .after=coup_id) %>%
  group_by(cowcode) %>%
  mutate(p_dem=case_when(
    is.na(y_post_coup) & (lag(e_polity2,n=1)>6 |
                            lag(e_polity2,n=2)>6 |
                            lag(e_polity2,n=3)>6 |
                            lag(e_polity2,n=4)>6)~1,
    TRUE~0
  )) %>%
  mutate(prior_dem= 
           (lag(e_polity2,n=5)+
              lag(e_polity2,n=4)+
              lag(e_polity2,n=3)+
              lag(e_polity2,n=2))/4) %>%
  ungroup %>%
  mutate(prior_dem=ifelse(prior_dem>=6,1,0))


## PWT controls ----

## Load PWT 10.0 from https://www.rug.nl/ggdc/productivity/pwt/pwt-releases/pwt100

pwt <- pwt %>%
  mutate(cowcode=countrycode(countrycode, origin = "iso3c", destination = "cown")) %>%
  filter(!is.na(cowcode)) %>%
  mutate(gdpc=rgdpe/pop) %>%
  select(cowcode, year, gdpc, pop)
data_main <- left_join(data_main, pwt, by=c("cowcode","year"))
data_1 <- left_join(data_1, pwt, by=c("cowcode","year"))
data_2 <- left_join(data_2, pwt, by=c("cowcode","year"))
data_4 <- left_join(data_4, pwt, by=c("cowcode","year"))
data_5 <- left_join(data_5, pwt, by=c("cowcode","year"))

## Constructing controls ----
data_main <- data_main %>%
  group_by(cowcode) %>%
  mutate(np_coup=ifelse(any_coup==1 & p_coup==0,1,0)) %>%
  mutate(l_torture=lag(torture,n=4)) %>%
  mutate(l_killings=lag(killings, n=4)) %>%
  mutate(l_repression=lag(repression, n=4)) %>%
  mutate(l_protest=lag(protest, n=4)) %>%
  mutate(cold_war=ifelse(year<1990,1,0)) %>%
  group_by(cowcode) %>%
  mutate(prior_coups=cumsum(coup)) %>%
  mutate(gdp_gr=(gdpc-lag(gdpc))/lag(gdpc)) %>%
  mutate(prior_repression=cumsum(replace_na(repression, 0))/row_number()) %>%
  mutate(prior_torture=cumsum(replace_na(torture, 0))/row_number()) %>%
  mutate(prior_killings=cumsum(replace_na(killings, 0))/row_number()) %>%
  ungroup

data_1 <- data_1 %>%
  group_by(cowcode) %>%
  mutate(np_coup=ifelse(any_coup==1 & p_coup==0,1,0)) %>%
  mutate(l_torture=lag(torture,n=4)) %>%
  mutate(l_killings=lag(killings, n=4)) %>%
  mutate(l_repression=lag(repression, n=4)) %>%
  mutate(l_protest=lag(protest, n=4)) %>%
  mutate(cold_war=ifelse(year<1990,1,0)) %>%
  group_by(cowcode) %>%
  mutate(prior_coups=cumsum(coup)) %>%
  mutate(gdp_gr=(gdpc-lag(gdpc))/lag(gdpc)) %>%
  mutate(prior_repression=cumsum(replace_na(repression, 0))/row_number()) %>%
  mutate(prior_torture=cumsum(replace_na(torture, 0))/row_number()) %>%
  mutate(prior_killings=cumsum(replace_na(killings, 0))/row_number()) %>%
  ungroup

data_2 <- data_2 %>%
  group_by(cowcode) %>%
  mutate(np_coup=ifelse(any_coup==1 & p_coup==0,1,0)) %>%
  mutate(l_torture=lag(torture,n=4)) %>%
  mutate(l_killings=lag(killings, n=4)) %>%
  mutate(l_repression=lag(repression, n=4)) %>%
  mutate(l_protest=lag(protest, n=4)) %>%
  mutate(cold_war=ifelse(year<1990,1,0)) %>%
  group_by(cowcode) %>%
  mutate(prior_coups=cumsum(coup)) %>%
  mutate(gdp_gr=(gdpc-lag(gdpc))/lag(gdpc)) %>%
  mutate(prior_repression=cumsum(replace_na(repression, 0))/row_number()) %>%
  mutate(prior_torture=cumsum(replace_na(torture, 0))/row_number()) %>%
  mutate(prior_killings=cumsum(replace_na(killings, 0))/row_number()) %>%
  ungroup

data_4 <- data_4 %>%
  group_by(cowcode) %>%
  mutate(np_coup=ifelse(any_coup==1 & p_coup==0,1,0)) %>%
  mutate(l_torture=lag(torture,n=4)) %>%
  mutate(l_killings=lag(killings, n=4)) %>%
  mutate(l_repression=lag(repression, n=4)) %>%
  mutate(l_protest=lag(protest, n=4)) %>%
  mutate(cold_war=ifelse(year<1990,1,0)) %>%
  group_by(cowcode) %>%
  mutate(prior_coups=cumsum(coup)) %>%
  mutate(gdp_gr=(gdpc-lag(gdpc))/lag(gdpc)) %>%
  mutate(prior_repression=cumsum(replace_na(repression, 0))/row_number()) %>%
  mutate(prior_torture=cumsum(replace_na(torture, 0))/row_number()) %>%
  mutate(prior_killings=cumsum(replace_na(killings, 0))/row_number()) %>%
  ungroup

data_5 <- data_5 %>%
  group_by(cowcode) %>%
  mutate(np_coup=ifelse(any_coup==1 & p_coup==0,1,0)) %>%
  mutate(l_torture=lag(torture,n=4)) %>%
  mutate(l_killings=lag(killings, n=4)) %>%
  mutate(l_repression=lag(repression, n=4)) %>%
  mutate(l_protest=lag(protest, n=4)) %>%
  mutate(cold_war=ifelse(year<1990,1,0)) %>%
  group_by(cowcode) %>%
  mutate(prior_coups=cumsum(coup)) %>%
  mutate(gdp_gr=(gdpc-lag(gdpc))/lag(gdpc)) %>%
  mutate(prior_repression=cumsum(replace_na(repression, 0))/row_number()) %>%
  mutate(prior_torture=cumsum(replace_na(torture, 0))/row_number()) %>%
  mutate(prior_killings=cumsum(replace_na(killings, 0))/row_number()) %>%
  ungroup


# Models ----
## Main ----
### Torture ----

m1_base <- lm(torture~
                p_coup +
                np_coup +
                l_torture, 
              data = data_main[!is.na(data_main$cowcode),])
cov1_b <- vcovCL(m1_base, cluster = ~cowcode)
r_se1_b <- sqrt(diag(cov1_b))

linearHypothesis(m1_base, "p_coup = np_coup", vcov. = cov1_b)

m1_base_c <- lm(torture~
           p_coup +
           np_coup +
           l_torture +
           coup_agency +
           casualties +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop,
         data=data_main)
cov1_bc <- vcovCL(m1_base_c, cluster = ~cowcode)
r_se1_bc <- sqrt(diag(cov1_bc))
summary_with_se <- coeftest(m1_base_c, vcov. = cov1_bc)
print(summary_with_se)

linearHypothesis(m1_base_c, "p_coup = np_coup", vcov. = cov1_bc)

m1 <- lm(torture~
           p_coup +
           np_coup +
           l_torture +
           coup_agency +
           casualties +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           protest +
           gdpc +
           gdp_gr,
         data=data_main)
cov1 <- vcovCL(m1, cluster = ~cowcode)
r_se1 <- sqrt(diag(cov1))

linearHypothesis(m1, "p_coup = np_coup", vcov. = cov1)

m2 <- lm(torture~
           any_coup +
           l_torture +
           coup_agency +
           casualties +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           protest +
           gdpc +
           gdp_gr,
         data=data_main)
cov2 <- vcovCL(m2, cluster = ~cowcode)
r_se2 <- sqrt(diag(cov2))


### Killings ----

m3_base <- lm(killings~
                p_coup +
                np_coup +
                l_killings, 
              data = data_main[!is.na(data_main$cowcode),])
cov3_b <- vcovCL(m3_base, cluster = ~cowcode)
r_se3_b <- sqrt(diag(cov3_b))

linearHypothesis(m3_base, "p_coup = np_coup", vcov. = cov3_b)

m3_base_c <- lm(killings~
                  p_coup +
                  np_coup +
                  l_killings +
                  coup_agency +
                  casualties +
                  e_polity2 +
                  p_dem +
                  cold_war +
                  prior_coups +
                  pop,
                data=data_main)
cov3_bc <- vcovCL(m3_base_c, cluster = ~cowcode)
r_se3_bc <- sqrt(diag(cov3_bc))
summary_with_se <- coeftest(m3_base_c, vcov. = cov3_bc)
print(summary_with_se)

linearHypothesis(m3_base_c, "p_coup = np_coup", vcov. = cov3_bc)

m3 <- lm(killings~
           p_coup +
           np_coup +
           l_killings +
           coup_agency +
           casualties +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           protest +
           gdpc +
           gdp_gr,
         data=data_main)
cov3 <- vcovCL(m3, cluster = ~cowcode)
r_se3 <- sqrt(diag(cov3))

linearHypothesis(m3, "p_coup = np_coup", vcov. = cov3)

m4 <- lm(killings~
           any_coup +
           l_killings +
           coup_agency +
           casualties +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           protest +
           gdpc +
           gdp_gr,
         data=data_main)
cov4 <- vcovCL(m4, cluster = ~cowcode)
r_se4 <- sqrt(diag(cov4))

### Repression ----
m5_base <- lm(repression~
                p_coup +
                np_coup +
                l_repression, 
              data = data_main[!is.na(data_main$cowcode),])
cov5_b <- vcovCL(m5_base, cluster = ~cowcode)
r_se5_b <- sqrt(diag(cov5_b))

linearHypothesis(m5_base, "p_coup = np_coup", vcov. = cov5_b)

m5_base_c <- lm(repression~
                  p_coup +
                  np_coup +
                  l_repression +
                  coup_agency +
                  casualties +
                  e_polity2 +
                  p_dem +
                  cold_war +
                  prior_coups +
                  pop,
                data=data_main)
cov5_bc <- vcovCL(m5_base_c, cluster = ~cowcode)
r_se5_bc <- sqrt(diag(cov5_bc))
summary_with_se <- coeftest(m5_base_c, vcov. = cov5_bc)
print(summary_with_se)

linearHypothesis(m5_base_c, "p_coup = np_coup", vcov. = cov5_bc)

m5 <- lm(repression~
           p_coup +
           np_coup +
           l_repression +
           coup_agency +
           casualties +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           protest +
           gdpc +
           gdp_gr,
         data=data_main)
cov5 <- vcovCL(m5, cluster = ~cowcode)
r_se5 <- sqrt(diag(cov5))

linearHypothesis(m5, "p_coup = np_coup", vcov. = cov5)

m6 <- lm(repression~
           any_coup +
           l_repression +
           coup_agency +
           casualties +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           protest +
           gdpc +
           gdp_gr,
         data=data_main)
cov6 <- vcovCL(m6, cluster = ~cowcode)
r_se6 <- sqrt(diag(cov6))

### Tables ----
stargazer(m1_base, m1_base_c, m1, m2,
          keep = c("p_coup",
                   "np_coup",
                   "any_coup",
                   "l_torture",
                   "coup_agency",
                   "casualties",
                   "e_polity2",
                   "p_dem",
                   "cold_war",
                   "prior_coups",
                   "pop",
                   "protest",
                   "gdpc",
                   "gdp_gr"),
          covariate.labels = c("Programmatic coup",
                               "Non-programmatic coup",
                               "Any coup",
                               "Prior torture",
                               "Coup agency",
                               "Casualties",
                               "Polity",
                               "Prior democracy",
                               "Cold War",
                               "Past coups",
                               "Population",
                               "Protest",
                               "GDP/capita",
                               "GDP growth"),
          se  = list(r_se1_b, r_se1_bc, r_se1, r_se2),
          keep.stat = c("n","rsq","F","adj.rsq"),
          table.layout = "-#-to-sa-n",
          title = "Coup types and torture",
          notes = "Robust standard errors clustered in country",
          header = F,
          df =F, 
          no.space = T,
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          type = "html",
          out = "table1.htm")

stargazer(m3_base, m3_base_c, m3,m4,
          keep = c("p_coup",
                   "np_coup",
                   "any_coup",
                   "l_killings",
                   "coup_agency",
                   "casualties",
                   "e_polity2",
                   "p_dem",
                   "cold_war",
                   "prior_coups",
                   "pop",
                   "protest",
                   "gdpc",
                   "gdp_gr"),
          covariate.labels = c("Programmatic coup",
                               "Non-programmatic coup",
                               "Any coup",
                               "Prior killings",
                               "Coup agency",
                               "Casualties",
                               "Polity",
                               "Prior democracy",
                               "Cold War",
                               "Past coups",
                               "Population",
                               "Protest",
                               "GDP/capita",
                               "GDP growth"),
          table.layout = "-#-t-sa-n",
          se  = list(r_se3_b, r_se3_bc, r_se3, r_se4),
          keep.stat = c("n","rsq","adj.rsq","F"),
          title = "Coup types and political killings",
          notes = "Robust standard errors clustered in country",
          header = F,
          df=F,
          no.space = T,
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          type = "html",
          out = "table2.htm")

stargazer(m5_base, m5_base_c, m5,m6,
          keep = c("p_coup",
                   "np_coup",
                   "any_coup",
                   "l_rep",
                   "coup_agency",
                   "casualties",
                   "e_polity2",
                   "p_dem",
                   "cold_war",
                   "prior_coups",
                   "pop",
                   "protest",
                   "gdpc",
                   "gdp_gr"),
          covariate.labels = c("Programmatic coup",
                               "Non-programmatic coup",
                               "Any coup",
                               "Past Repression",
                               "Coup agency",
                               "Casualties",
                               "Polity",
                               "Prior democracy",
                               "Cold War",
                               "Past coups",
                               "Population",
                               "Protest",
                               "GDP/capita",
                               "GDP growth"),
          table.layout = "-#-t-sa-n",
          se = list(r_se5_b, r_se5_bc, r_se5, r_se6),
          keep.stat = c("n","rsq","adj.rsq","F"),
          title = "Coup types and repression",
          notes = "Robust standard errors clustered in country",
          header = F,
          df=F,
          no.space = T,
          digits = 2,
          digits.extra = 2,
          type = "html",
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          out = "table3.htm")

#### Interpretation of effects in text ----
m1$coefficients[2]/sd(data_main$torture)
m3$coefficients[2]/sd(data_main$killings)
m5$coefficients[2]/sd(data_main$repression, na.rm = T)

m1$coefficients[3]/sd(data_main$torture)
m3$coefficients[3]/sd(data_main$killings)
m5$coefficients[3]/sd(data_main$repression, na.rm = T)

#### Coefficient plot ----

coef <- data.frame("iv"=c("Programmatic","Non-programmatic","Any",
                          "Programmatic","Non-programmatic","Any",
                          "Programmatic","Non-programmatic","Any"),
                   "dv"=c("Torture","Torture","Torture",
                          "Killings","Killings","Killings",
                          "Repression","Repression","Repression"),
                   "coef"=rep(NA,9),
                   "se"=rep(NA,9))

coef$coef[1] <- m1$coefficients[2]
coef$coef[2] <- m1$coefficients[3]
coef$coef[3] <- m2$coefficients[2]
coef$coef[4] <- m3$coefficients[2]
coef$coef[5] <- m3$coefficients[3]
coef$coef[6] <- m4$coefficients[2]
coef$coef[7] <- m5$coefficients[2]
coef$coef[8] <- m5$coefficients[3]
coef$coef[9] <- m6$coefficients[2]
coef$se[1] <- r_se1[2]
coef$se[2] <- r_se1[3]
coef$se[3] <- r_se2[2]
coef$se[4] <- r_se3[2]
coef$se[5] <- r_se3[3]
coef$se[6] <- r_se4[2]
coef$se[7] <- r_se5[2]
coef$se[8] <- r_se5[3]
coef$se[9] <- r_se6[2]
coef$lower95 <- coef$coef-(1.96*coef$se)
coef$upper95 <- coef$coef+(1.96*coef$se)
coef$lower90 <- coef$coef-(1.645*coef$se)
coef$upper90 <- coef$coef+(1.645*coef$se)

coef <- coef %>%
  mutate(iv=factor(iv, levels = c("Any",
                                  "Non-programmatic",
                                  "Programmatic"))) %>%
  mutate(dv=factor(dv, levels = c("Torture",
                                  "Killings",
                                  "Repression"))) %>%
  mutate(order=case_when(
    dv=="Torture"~1,
    dv=="Killings"~2,
    dv=="Repression"~3
  ))

line_styles <- c("dashed", "dotted", "solid")  # Choose the line styles you want
point_symbols <- c(16, 17, 18)  # Choose the point symbols you want

plot <- ggplot(coef, aes(y=coef, x=iv, color=reorder(dv,-order))) +
  geom_point(position=position_dodge(0.9), 
             aes(shape=reorder(dv,-order)), 
             size=3) +
  geom_errorbar(aes(ymin=lower95, 
                    ymax=upper95),
                width=0,
                position=position_dodge(0.9)) +
  geom_errorbar(aes(ymin=lower90, 
                    ymax=upper90),
                width=0,
                size=1,
                alpha=0.6,
                position=position_dodge(0.9)) +
  geom_hline(yintercept=0, linetype="dashed") +
  theme_classic(base_family = "Times New Roman") +
  theme(text = element_text(size = 22,
                            color = "black"),
        axis.text = element_text(color = "black"),  
        axis.title = element_text(color = "black")) +
  coord_flip() +
  xlab("") +
  ylab("") +
  scale_color_manual(values = c("black","black","black"), 
                       guide = guide_legend(title="", reverse = T)) +
  scale_shape_manual(values = point_symbols,
                     guide = guide_legend(title="", reverse = T))  
plot
ggsave("coefplot_lines.jpg", 
       dpi=300,
       width=9,
       height=6)

### Protest ----
m1_base <- lm(protest~
           p_coup +
           np_coup +
           l_protest,
         data=data_main[!is.na(data_main$cowcode),])
cov1_base <- vcovCL(m1_base, cluster = ~cowcode)
r_se1_b <- sqrt(diag(cov1_base))

m1_base_c <- lm(protest~
           p_coup +
           np_coup +
           l_protest +
           coup_agency +
           casualties +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop,
         data=data_main)
cov1_base_c <- vcovCL(m1_base_c, cluster = ~cowcode)
r_se1_bc <- sqrt(diag(cov1_base_c))

m1 <- lm(protest~
            p_coup +
            np_coup +
            l_protest +
            coup_agency +
            casualties +
            e_polity2 +
            p_dem +
            cold_war +
            prior_coups +
            pop +
            gdpc +
            gdp_gr,
          data=data_main)
cov1 <- vcovCL(m1, cluster = ~cowcode)
r_se1 <- sqrt(diag(cov1))

m2 <- lm(protest~
            any_coup +
            l_protest +
            coup_agency +
            casualties +
            e_polity2 +
            p_dem +
            cold_war +
            prior_coups +
            pop +
            gdpc +
            gdp_gr,
          data=data_main)
cov2 <- vcovCL(m2, cluster = ~cowcode)
r_se2 <- sqrt(diag(cov2))

#### Table ----
stargazer(m1_base, m1_base_c, m1, m2,
          keep = c("p_coup",
                   "np_coup",
                   "any_coup",
                   "l_protest",
                   "coup_agency",
                   "casualties",
                   "e_polity2",
                   "p_dem",
                   "cold_war",
                   "prior_coups",
                   "pop",
                   "gdpc",
                   "gdp_gr"),
          covariate.labels = c("Programmatic coup",
                               "Non-programmatic coup",
                               "Any coup",
                               "Prior protest",
                               "Coup agency",
                               "Casualties",
                               "Polity",
                               "Prior democracy",
                               "Cold War",
                               "Past coups",
                               "Population",
                               "GDP/capita",
                               "GDP growth"),
          se  = list(r_se1_b, r_se1_bc, r_se1, r_se2),
          keep.stat = c("n","rsq","adj.rsq","F"),
          table.layout = "-#-to-sa-n",
          title = "Coup types and protest",
          notes = "Robust standard errors clustered in country",
          header = F,
          df =F, 
          no.space = T,
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          type = "html",
          out = "table_protest.htm")

# Appendix ----
## 1-year window ----
### Torture ----
t1y1 <- lm(torture~
           p_coup +
           np_coup +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_1)
cov1 <- vcovCL(t1y1, cluster = ~cowcode)
se_t1y1 <- sqrt(diag(cov1))

t2y1 <- lm(torture~
           any_coup +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_1)
cov2 <- vcovCL(t2y1, cluster = ~cowcode)
se_t2y1 <- sqrt(diag(cov2))

### Killings ----
k1y1 <- lm(killings~
           p_coup +
           np_coup +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_1)
cov3 <- vcovCL(k1y1, cluster = ~cowcode)
se_k1y1 <- sqrt(diag(cov3))

k2y1 <- lm(killings~
           any_coup +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_1)
cov4 <- vcovCL(k2y1, cluster = ~cowcode)
se_k2y1 <- sqrt(diag(cov4))

### Repression ----
r1y1 <- lm(repression~
           p_coup +
           np_coup +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_1)
cov5 <- vcovCL(r1y1, cluster = ~cowcode)
se_r1y1 <- sqrt(diag(cov5))

r2y1 <- lm(repression~
           any_coup +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_1)
cov6 <- vcovCL(r2y1, cluster = ~cowcode)
se_r2y1 <- sqrt(diag(cov6))

## 2-year window ----
### Torture ----
t1y2 <- lm(torture~
           p_coup +
           np_coup +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_2)
cov1 <- vcovCL(t1y2, cluster = ~cowcode)
se_t1y2 <- sqrt(diag(cov1))


t2y2 <- lm(torture~
           any_coup +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_2)
cov2 <- vcovCL(t2y2, cluster = ~cowcode)
se_t2y2 <- sqrt(diag(cov2))

### Killings ----
k1y2 <- lm(killings~
           p_coup +
           np_coup +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_2)
cov3 <- vcovCL(k1y2, cluster = ~cowcode)
se_k1y2 <- sqrt(diag(cov3))

k2y2 <- lm(killings~
           any_coup +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_2)
cov4 <- vcovCL(k2y2, cluster = ~cowcode)
se_k2y2 <- sqrt(diag(cov4))

### Repression ----
r1y2 <- lm(repression~
           p_coup +
           np_coup +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_2)
cov5 <- vcovCL(r1y2, cluster = ~cowcode)
se_r1y2 <- sqrt(diag(cov5))

r2y2 <- lm(repression~
           any_coup +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_2)
cov6 <- vcovCL(r2y2, cluster = ~cowcode)
se_r2y2 <- sqrt(diag(cov6))

## 4-year window ----
### Torture ----
t1y4 <- lm(torture~
           p_coup +
           np_coup +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_4)
cov1 <- vcovCL(t1y4, cluster = ~cowcode)
se_t1y4 <- sqrt(diag(cov1))


t2y4 <- lm(torture~
           any_coup +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_4)
cov2 <- vcovCL(t2y4, cluster = ~cowcode)
se_t2y4 <- sqrt(diag(cov2))

### Killings ----
k1y4 <- lm(killings~
           p_coup +
           np_coup +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_4)
cov3 <- vcovCL(k1y4, cluster = ~cowcode)
se_k1y4 <- sqrt(diag(cov3))

k2y4 <- lm(killings~
           any_coup +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_4)
cov4 <- vcovCL(k2y4, cluster = ~cowcode)
se_k2y4 <- sqrt(diag(cov4))

### Repression ----
r1y4 <- lm(repression~
           p_coup +
           np_coup +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_4)
cov5 <- vcovCL(r1y4, cluster = ~cowcode)
se_r1y4 <- sqrt(diag(cov5))

r2y4 <- lm(repression~
           any_coup +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_4)
cov6 <- vcovCL(r2y4, cluster = ~cowcode)
se_r2y4 <- sqrt(diag(cov6))

## 5-year window ----
### Torture ----
t1y5 <- lm(torture~
           p_coup +
           np_coup +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_5)
cov1 <- vcovCL(t1y5, cluster = ~cowcode)
se_t1y5 <- sqrt(diag(cov1))


t2y5 <- lm(torture~
           any_coup +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_5)
cov2 <- vcovCL(t2y5, cluster = ~cowcode)
se_t2y5 <- sqrt(diag(cov2))

### Killings ----
k1y5 <- lm(killings~
           p_coup +
           np_coup +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_5)
cov3 <- vcovCL(k1y5, cluster = ~cowcode)
se_k1y5 <- sqrt(diag(cov3))

k2y5 <- lm(killings~
           any_coup +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_5)
cov4 <- vcovCL(k2y5, cluster = ~cowcode)
se_k2y5 <- sqrt(diag(cov4))

### Repression ----
r1y5 <- lm(repression~
           p_coup +
           np_coup +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_5)
cov5 <- vcovCL(r1y5, cluster = ~cowcode)
se_r1y5 <- sqrt(diag(cov5))

r2y5 <- lm(repression~
           any_coup +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_5)
cov6 <- vcovCL(r2y5, cluster = ~cowcode)
se_r2y5 <- sqrt(diag(cov6))

### Tables ----
stargazer(t1y1, t2y2, t1y2, t2y2, t1y4, t2y4, t1y5, t2y5,
          keep = c("p_coup",
                   "np_coup",
                   "any_coup",
                   "l_torture",
                   "coup_agency",
                   "casualties",
                   "protest",
                   "e_polity2",
                   "p_dem",
                   "cold_war",
                   "prior_coups",
                   "pop",
                   "gdpc",
                   "gdp_gr"),
          covariate.labels = c("Programmatic coup",
                               "Non-programmatic coup",
                               "Any coup",
                               "Prior torture",
                               "Coup agency",
                               "Casualties",
                               "Protest",
                               "Polity",
                               "Prior democracy",
                               "Cold War",
                               "Past coups",
                               "Population",
                               "GDP/capita",
                               "GDP growth"),
          se  = list(se_t1y1, se_t2y2, se_t1y2, se_t2y2, se_t1y4, se_t2y4, se_t1y5, se_t2y5),
          keep.stat = c("n","rsq","adj.rsq"),
          table.layout = "-#-to-sa-n",
          title = "Coup types and torture",
          notes = "Robust standard errors clustered in country",
          header = F,
          df =F, 
          no.space = T,
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          type = "html",
          out = "table1_other_windows.htm")

stargazer(k1y1, k2y2, k1y2, k2y2, k1y4, k2y4, k1y5, k2y5,
          keep = c("p_coup",
                   "np_coup",
                   "any_coup",
                   "l_killings",
                   "coup_agency",
                   "casualties",
                   "protest",
                   "e_polity2",
                   "p_dem",
                   "cold_war",
                   "prior_coups",
                   "pop",
                   "gdpc",
                   "gdp_gr"),
          covariate.labels = c("Programmatic coup",
                               "Non-programmatic coup",
                               "Any coup",
                               "Prior killings",
                               "Coup agency",
                               "Casualties",
                               "Protest",
                               "Polity",
                               "Prior democracy",
                               "Cold War",
                               "Past coups",
                               "Population",
                               "GDP/capita",
                               "GDP growth"),
          table.layout = "-#-t-sa-n",
          se  = list(se_k1y1, se_k2y2, se_k1y2, se_k2y2, se_k1y4, se_k2y4, se_k1y5, se_k2y5),
          keep.stat = c("n","rsq","adj.rsq"),
          title = "Coup types and political killings",
          notes = "Robust standard errors clustered in country",
          header = F,
          df=F,
          no.space = T,
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          type = "html",
          out = "table2_other_windows.htm")

stargazer(r1y1, r2y2, r1y2, r2y2, r1y4, r2y4, r1y5, r2y5,
          keep = c("p_coup",
                   "np_coup",
                   "any_coup",
                   "l_rep",
                   "coup_agency",
                   "casualties",
                   "protest",
                   "e_polity2",
                   "p_dem",
                   "cold_war",
                   "prior_coups",
                   "pop",
                   "gdpc",
                   "gdp_gr"),
          covariate.labels = c("Programmatic coup",
                               "Non-programmatic coup",
                               "Any coup",
                               "Past Repression",
                               "Coup agency",
                               "Casualties",
                               "Protest",
                               "Polity",
                               "Prior democracy",
                               "Cold War",
                               "Past coups",
                               "Population",
                               "GDP/capita",
                               "GDP growth"),
          table.layout = "-#-t-sa-n",
          se  = list(se_r1y1, se_r2y2, se_r1y2, se_r2y2, se_r1y4, se_r2y4, se_r1y5, se_r2y5),
          keep.stat = c("n","rsq","adj.rsq"),
          title = "Coup types and repression",
          notes = "Robust standard errors clustered in country",
          header = F,
          df=F,
          no.space = T,
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          type = "html",
          out = "table3_other_windows.htm")

## only post-coup years ----

m1_coups <- lm(torture~
                 p_coup +
                 l_torture, 
               data = data_main[!is.na(data_main$cowcode) & data_main$any_coup==1,])
cov1_c <- vcovCL(m1_coups, cluster = ~cowcode)
r_se1_c <- sqrt(diag(cov1_c))
summary_with_se <- coeftest(m1_coups, vcov. = cov1_c)
print(summary_with_se)

m2_coups <- lm(killings~
                 p_coup +
                 l_killings, 
               data = data_main[!is.na(data_main$cowcode) & data_main$any_coup==1,])
cov2_c <- vcovCL(m2_coups, cluster = ~cowcode)
r_se2_c <- sqrt(diag(cov2_c))
summary_with_se <- coeftest(m2_coups, vcov. = cov2_c)
print(summary_with_se)

m3_coups <- lm(repression~
                 p_coup +
                 l_repression, 
               data = data_main[!is.na(data_main$cowcode) & data_main$any_coup==1,])
cov3_c <- vcovCL(m3_coups, cluster = ~cowcode)
r_se3_c <- sqrt(diag(cov3_c))
summary_with_se <- coeftest(m3_coups, vcov. = cov3_c)
print(summary_with_se)

## Three-level factor ----

data_main <- data_main %>%
  mutate(coup_types=case_when(
    p_coup==1~"programmatic",
    np_coup==1~"non-programmatic",
    TRUE~"no coup"
  ),
  coup_types = fct_relevel(factor(coup_types), "non-programmatic")) 

m1 <- lm(torture~
           coup_types +
           l_torture +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_main)
cov1 <- vcovCL(m1, cluster = ~cowcode)
r_se1 <- sqrt(diag(cov1))
summary(m1)

### Killings ----
m2 <- lm(killings~
           coup_types +
           l_killings +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_main)
cov2 <- vcovCL(m2, cluster = ~cowcode)
r_se2 <- sqrt(diag(cov2))
summary_with_se <- coeftest(m2, vcov. = cov2)
print(summary_with_se)

### Repression ----
m3 <- lm(repression~
           coup_types +
           l_repression +
           coup_agency +
           casualties +
           protest +
           e_polity2 +
           p_dem +
           cold_war +
           prior_coups +
           pop +
           gdpc +
           gdp_gr,
         data=data_main)
cov3 <- vcovCL(m3, cluster = ~cowcode)
r_se3 <- sqrt(diag(cov3))

### Tables ----
stargazer(m1_coups,m2_coups,m3_coups,
          keep = c("p_coup",
                   "l_torture",
                   "l_killings",
                   "l_repression"),
          covariate.labels = c("Programmatic coup",
                               "Lagged DV",
                               "Lagged DV",
                               "Lagged DV"),
          keep.stat = c("n","rsq","F","adj.rsq"),
          table.layout = "-#-to-sa-n",
          se  = list(r_se1_c, r_se2_c, r_se3_c),
          title = "Only post-coup observations",
          notes = "Robust standard errors clustered in country",
          header = F,
          df =F, 
          no.space = T,
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          type = "html",
          out = "table_postcoup.htm")

stargazer(m1,m2,m3,
          keep = c("coup_types",
                   "l_torture",
                   "l_killings",
                   "l_repression",
                   "coup_agency",
                   "casualties",
                   "protest",
                   "e_polity2",
                   "p_dem",
                   "cold_war",
                   "prior_coups",
                   "pop",
                   "gdpc",
                   "gdp_gr"),
          covariate.labels = c("No coup",
                               "Programmatic coup",
                               "Lagged DV",
                               "Lagged DV",
                               "Lagged DV",
                               "Coup agency",
                               "Casualties",
                               "Protest",
                               "Polity",
                               "Prior democracy",
                               "Cold War",
                               "Past coups",
                               "Population",
                               "GDP/capita",
                               "GDP growth"),
          keep.stat = c("n","rsq","F","adj.rsq"),
          table.layout = "-#-to-sa-n",
          se  = list(r_se1, r_se2, r_se3),
          title = "Non-programmatic coups as excluded category",
          notes = "Robust standard errors clustered in country",
          header = F,
          df =F, 
          no.space = T,
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("†", "*", "**", "***"),
          type = "html",
          out = "table_factor.htm")