################################################################################
####Replication code for paper################################################
################################################################################

library(foreign)
library(apsrtable)
library(lmtest)
library(sandwich)
library(pcse)
library(car)
library(gee)

setwd("/Users/frj/Documents/Gammelt/Berkeley/PartyOrg_paper")
X<-read.csv(file="Datafiles/Data_paper.csv")

attach(X)

##Key variables
# enop_sea = Effective number of parties, calculated by seats held in the state assembly
# enop_vot = Effective number of parties, calculated by votes
# ev_seats = electoral volatility, calculated by seats
# ev_votes = electoral volatility, calculated by votes
# ev_seats_0 = electoral volatility by seats in election t-1
# ev_votes_0 = electoral volatility by votes in election t-1
# enop_vo0 = Effective number of parties, calculated by votes, in election t-1
# strparty = parties coded as organized, coded 3
# strparty23 = parties coded as organized, coded as 2 or 3
# proporg = proportion of parties that were coded as 3 organized
# proporg23 = proportion of parties in state coded as 2 or 3 organized
# ENgr2004 = Effective number of groups in the state, coded from NES 2004

#FIGURE 1
electionlevels<-levels(as.factor(electio0))
ENOP<-as.data.frame(matrix(nrow=length(electionlevels), ncol=3))
for (i in 1:length(electionlevels)) {
ENOP[i, 1]<-electionlevels[i]
ENOP[i, 2]<-mean(enop_sea[electio0==electionlevels[i]], na.rm=T)
ENOP[i, 3]<-mean(enop_vot[electio0==electionlevels[i]], na.rm=T)
	}
pdf(file="Fig1.pdf")
plot(ENOP$V2~ENOP$V1, type="l", col="blue", xlab="Year", ylab="Effective number of parties", ylim=c(1.5, 7))
names<-c("Eff. number of parties by seats", "Eff. number og parties by votes")
legend("topright", names, col=c("blue", "black"), lty=c(1, 8), cex=0.8)
par(new=TRUE)
plot(ENOP$V3~ENOP$V1, type="l", lty=8, col="black", ylim=c(1.5, 7), ylab="", xlab="")
dev.off()

#######FIGURE 2
statelevels<-levels(as.factor(state))
state_av<-as.data.frame(matrix(nrow=length(statelevels), ncol=6))

for (i in 1:length(statelevels)) {
state_av[i,1]<-statelevels[i]
state_av[i, 2]<-enop_vot[state==statelevels[i] & election==1]
state_av[i, 3]<-enop_vot[state==statelevels[i] & election==8]
state_av[i, 4]<-enop_sea[state==statelevels[i] & election==1]
state_av[i, 5]<-enop_sea[state==statelevels[i] & election==8]
}

names.arg<-ifelse(state_av$V1=="Andhra Pradesh", "AP", NA)
names.arg<-ifelse(state_av$V1=="Assam", "AS", names.arg)
names.arg<-ifelse(state_av$V1=="Bihar", "BI", names.arg)
names.arg<-ifelse(state_av$V1=="Gujarat", "GUJ", names.arg)
names.arg<-ifelse(state_av$V1=="Haryana", "HAR", names.arg)
names.arg<-ifelse(state_av$V1=="Karnataka", "KN", names.arg)
names.arg<-ifelse(state_av$V1=="Kerala", "KER", names.arg)
names.arg<-ifelse(state_av$V1=="Madhya Pradesh", "MP", names.arg)
names.arg<-ifelse(state_av$V1=="Maharashtra", "MAH", names.arg)
names.arg<-ifelse(state_av$V1=="Orissa", "OR", names.arg)
names.arg<-ifelse(state_av$V1=="Punjab", "PUN", names.arg)
names.arg<-ifelse(state_av$V1=="Rajasthan", "RAJ", names.arg)
names.arg<-ifelse(state_av$V1=="Tamil Nadu", "TN", names.arg)
names.arg<-ifelse(state_av$V1=="Uttar Pradesh", "UP", names.arg)
names.arg<-ifelse(state_av$V1=="West Bengal", "WB", names.arg)
names.arg<-names.arg[!is.na(names.arg)]

state_av[,6]<-names.arg

names(state_av)<-c("State", "ENOP_vot1", "ENOP_vot8", "ENOP_sea1", "ENOP_sea8", "state_short")

col<-c("gray", "brown")
names<-c("First election", "Eight election")

for (i in 1:length(statelevels)) {
state_av[i,1]<-statelevels[i]
state_av[i, 2]<-mean(enop_vot[state==statelevels[i]], na.rm=T)
state_av[i, 3]<-sd(enop_vot[state==statelevels[i]], na.rm=T)
state_av[i, 4]<-mean(enop_sea[state==statelevels[i]], na.rm=T)
state_av[i, 5]<-sd(enop_sea[state==statelevels[i]], na.rm=T)
}

pdf(file="Fig2.pdf")
par(oma=c(2.5,0.8,0.5,0.5))
boxplot(enop_vot~state, names.arg=statelevels, cex.names=0.8, beside=T, col=c("lightgrey"),
ylab="Effective number of parties, by votes",
xlab="", 
ylim=c(2,9),
cex.axis=0.8,
cex.main=1, las=2
)
dev.off()

######Figure 3

pdf(file="Fig3.pdf")
par(mfrow=c(2,2))
plot(as.factor(nparties), main="Number of parties", ylab="Frequency", ylim=c(0, 60))
plot(density(enop_vot), main="Effective number of parties", ylab="Density", xlab="")
plot(factor(strparty, levels=c(0,1,2,3,4,5,6, 7)), main="Number of organized parties", ylab="Frequency", ylim=c(0, 80))
plot(density(proporg), xlim=c(0,1), main="Proportion organized parties", ylab="Density", xlab="")
dev.off()

####Figure 4

pdf(file="Fig4.pdf", height=4)
par(mfrow=c(1,2))
plot(enop_vot~proporg, xlab="Prop. organized parties", ylab="Eff. no of parties", main="Effective number of parties")
lines(lowess(enop_vot~proporg), col = 4, lty=3)
legend("topright", "Loess curve", col = 4, lty=3)

plot(ev_votes[complete.cases(ev_votes)]~proporg[complete.cases(ev_votes)], xlab="Prop. organized parties", ylab="Electoral volatility", main="Electoral volatility")
lines(lowess(ev_votes[complete.cases(ev_votes)]~proporg[complete.cases(ev_votes)], f=1.5), col = 4, lty=3)
legend("topright", "Loess curve", col = 4, lty=3)
dev.off()

#####Figure 5
MP<-c(64.38, 76.64, 97.8, 88.2, 90.9, 50)
UP<-c(5.3, 13.4, 75.3, 0, 10.2, 31.6)
pdf(file="Fig5.pdf", height=4)
par(mfrow=c(1,2))
par(oma=c(1.5,0.8,0.5,0.5))
barplot(MP, names.arg=c("Upper Caste", "OBCs", "SCs", "STs", "Muslims", "Others"), 
main="Support for INC in Madhya Pradesh 1980", ylim=c(0,100), cex.main=0.8, las=2)
barplot(UP, names.arg=c("Upper Caste", "OBCs", "SCs", "STs", "Muslims", "Others"), 
main="Support for BSP in Uttar Pradesh in 2004 ", ylim=c(0,75), cex.main=0.8, las=2)
dev.off()

###################
#######OLS models for ENOP

model1<-lm(enop_vot~enop_vo0 +proporg+I(proporg^2)+ev_votes, data=X)

model2<-lm(enop_vot~enop_vo0 +proporg+I(proporg^2)+ev_votes+ENgr2004+socialba+lnrealin+lnfiscal+turnout+yr1991,data=X)

model3<-lm(enop_vot~enop_vo0 +proporg+I(proporg^2)+ev_votes+socialba+cleavage+lnrealin+lnfiscal+turnout+yr1991, data=X)

model4<-lm(enop_vot~enop_vo0 +proporg+I(proporg^2)+ev_votes+socialba+eff_clus+lnrealin+lnfiscal+turnout+yr1991, data=X)

model5<-lm(enop_vot~proporg+I(proporg^2)+ev_votes, data=X)

model6<-lm(enop_vot~proporg+I(proporg^2)+ev_votes+ENgr2004+socialba+lnrealin+lnfiscal+turnout+yr1991,data=X)

model7<-lm(enop_vot~proporg+I(proporg^2)+ev_votes+socialba+cleavage+lnrealin+lnfiscal+turnout+yr1991, data=X)

model8<-lm(enop_vot~proporg+I(proporg^2)+ev_votes+socialba+eff_clus+lnrealin+lnfiscal+turnout+yr1991, data=X)

vif(model2)
vif(model3)
vif(model4)

##Panel corrected standard errors

model1$se<-vcovPC(model1, groupN=id[!is.na(ev_votes)], groupT=election[!is.na(ev_votes)])
model2$se<-vcovPC(model2, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)])
model3$se<-vcovPC(model3, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, cleavage, lnrealin, lnfiscal, turnout, yr1991)])

model4$se<-vcovPC(model4, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, eff_clus, lnrealin, lnfiscal, turnout, yr1991)])

model5$se<-vcovPC(model5, groupN=id[!is.na(ev_votes)], groupT=election[!is.na(ev_votes)])
model6$se<-vcovPC(model6, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)])
model7$se<-vcovPC(model7, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, cleavage, lnrealin, lnfiscal, turnout, yr1991)])

model8$se<-vcovPC(model8, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, eff_clus, lnrealin, lnfiscal, turnout, yr1991)])

table_OLS<-apsrtable(model1, model2, model3, model4, digits=2, stars="default", coef.names=c("Intercept", "Lagged eff. number of parties", "Proportion organized parties", "Proportion organized parties squared", "Electoral volatility", "Effective number of groups 2004", "Dummy for social base", "Log (Real Income)", "Log (Fiscal spending)", "Electoral turnout", "Dummy for post 1991", "Cleavages", "Effective number of clusters"), se="robust")

table_nolag<-apsrtable(model5, model6, model7, model8, digits=2, stars="default", coef.names=c("Intercept", "Proportion organized parties", "Proportion organized parties squared", "Electoral volatility", "Effective number of groups 2004", "Dummy for social base", "Log (Real Income)", "Log (Fiscal spending)", "Electoral turnout", "Dummy for post 1991", "Cleavages", "Effective number of clusters"), se="robust")

##Output Table 1
table_OLS

#####GEE models

model9<-gee(enop_vot~enop_vo0 +proporg+I(proporg^2)+ev_votes, id=id, data=X, corstr="exchangeable", maxiter=100)

model10<-gee(enop_vot~enop_vo0 +proporg+I(proporg^2)+ev_votes+ENgr2004+socialba+lnrealin+lnfiscal+turnout+yr1991, id=id, data=X, corstr="exchangeable", maxiter=100)

model11<-gee(enop_vot~enop_vo0 +proporg+I(proporg^2)+ev_votes+socialba+cleavage+lnrealin+lnfiscal+turnout+yr1991, id=id, data=X, corstr="exchangeable", maxiter=100)

model12<-gee(enop_vot~enop_vo0 +proporg+I(proporg^2)+ev_votes+socialba+eff_clus+lnrealin+lnfiscal+turnout+yr1991, id=id, data=X, corstr="exchangeable", maxiter=100)

table_GEE<-apsrtable(model9, model10, model11, model12, digits=2, stars="default", coef.names=c("Intercept", "Lagged eff. number of parties", "Proportion organized parties", "Prop. organized parties squared", "Electoral volatility", "Eff. number of groups 2004", "Dummy for social base", "Log (Real Income)", "Log (Fiscal spending)", "Electoral turnout", "Dummy for post 1991", "Cleavages", "Effective number of clusters"))

##Output Appendix A
table_GEE

#######Electoral volatility

model1<-lm(ev_votes~ev_votes_0+proporg, data=X)

model2<-lm(ev_votes~ev_votes_0+proporg+enop_vot+ENgr2004+socialba+lnrealin+lnfiscal+turnout+yr1991, data=X)

model3<-lm(ev_votes~ev_votes_0+proporg+socialba+cleavage+lnrealin+lnfiscal+turnout+yr1991, data=X)

model4<-lm(ev_votes~ev_votes_0+proporg+socialba+eff_clus+lnrealin+lnfiscal+turnout+yr1991, data=X)

model5<-lm(ev_votes~proporg, data=X)

model6<-lm(ev_votes~proporg+enop_vot+ENgr2004+socialba+lnrealin+lnfiscal+turnout+yr1991,data=X)

model7<-lm(ev_votes~proporg+socialba+cleavage+lnrealin+lnfiscal+turnout+yr1991, data=X)

model8<-lm(ev_votes~proporg+socialba+eff_clus+lnrealin+lnfiscal+turnout+yr1991, data=X)

vif(model2)
vif(model3)
vif(model4)

##Panel corrected standard errors
model1$se<-vcovPC(model1, groupN=id[!is.na(ev_votes_0)], groupT=election[!is.na(ev_votes_0)])
model2$se<-vcovPC(model2, groupN=id[complete.cases(ev_votes_0, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes_0, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)])
model3$se<-vcovPC(model3, groupN=id[complete.cases(ev_votes_0, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes_0, ENgr2004, cleavage, lnrealin, lnfiscal, turnout, yr1991)])

model4$se<-vcovPC(model4, groupN=id[complete.cases(ev_votes_0, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes_0, ENgr2004, eff_clus, lnrealin, lnfiscal, turnout, yr1991)])

model5$se<-vcovPC(model5, groupN=id[!is.na(ev_votes)], groupT=election[!is.na(ev_votes)])
model6$se<-vcovPC(model6, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)])
model7$se<-vcovPC(model7, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, cleavage, lnrealin, lnfiscal, turnout, yr1991)])

model8$se<-vcovPC(model8, groupN=id[complete.cases(ev_votes, ENgr2004, socialba, lnrealin, lnfiscal, turnout, yr1991)], groupT=election[complete.cases(ev_votes, ENgr2004, eff_clus, lnrealin, lnfiscal, turnout, yr1991)])

table_EV<-apsrtable(model1, model2, model3, model4, digits=2, stars="default", coef.names=c("Intercept", "Lagged e. volatility", "Proportion organized parties", "Effective number of parties", "Effective number of groups 2004", "Dummy for social base", "Log (Real Income)", "Log (Fiscal spending)", "Electoral turnout", "Dummy for post 1991", "Cleavages", "Effective number of clusters"), se="robust")

table_EV_nolag<-apsrtable(model5, model6, model7, model8, digits=2, stars="default", coef.names=c("Intercept","Proportion organized parties", "Effective number of parties", "Effective number of groups 2004", "Dummy for social base", "Log (Real Income)", "Log (Fiscal spending)", "Electoral turnout", "Dummy for post 1991", "Cleavages", "Effective number of clusters"), se="robust")

##Output Table 2
table_EV

#####VOLATILITY GEE MODELS

model9<-gee(ev_votes~ev_votes_0+proporg, id=id, data=X, corstr="exchangeable", maxiter=100)

model10<-gee(ev_votes~ev_votes_0+proporg+enop_vot+ENgr2004+socialba+lnrealin+lnfiscal+turnout+yr1991,id=id, data=X, corstr="exchangeable", maxiter=100)

model11<-gee(ev_votes~ev_votes_0+proporg+socialba+cleavage+lnrealin+lnfiscal+turnout+yr1991, id=id, data=X, corstr="exchangeable", maxiter=100)

model12<-gee(ev_votes~ev_votes_0+proporg+socialba+eff_clus+lnrealin+lnfiscal+turnout+yr1991, id=id, data=X, corstr="exchangeable", maxiter=100)

table_EV_GEE<-apsrtable(model9, model10, model11, model12, digits=2, stars="default", coef.names=c("Intercept", "Lagged e. volatility", "Proportion organized parties", "Effective number of parties", "Effective number of groups 2004", "Dummy for social base", "Log (Real Income)", "Log (Fiscal spending)", "Electoral turnout", "Dummy for post 1991", "Cleavages", "Effective number of clusters"))

## Appendix table 3
table_EV_GEE