## ## Drew A. Linzer ## dlinzer@emory.edu ## http://userwww.service.emory.edu/~dlinzer ## ## May 27, 2007 ## ## LR-world_prices_analysis.R ## ## R code to replicate the statistical analyses in: ## ## Drew A. Linzer and Ronald L. Rogowski. 2008. "Lower Prices: ## The Impact of Majoritarian Systems in Democracies Around the ## World." The Journal of Politics. 70(1): 17-27. ## ## http://journals.cambridge.org/action/displayIssue?jid=JOP&volumeId=70&issueId=01 ## ## This replication file contains commands to run analyses that ## appear in both the published paper and in the web appendix. ## library(gee) library(Amelia) library(plm) options(max.print=999999) par(mar=c(5,4.5,0.5,0.5),mfrow=c(1,1),ask=F) setwd("C:/.../") # Enter working directory here ## ## Load in data file: select either Freedom House (FH) or Polity. ## whichdata <- "FH" #whichdata <- "Polity" fullmerge <- read.csv("Full-merged-data.csv") if (whichdata=="FH") { cat("\n ##############################################################") cat("\n ## USING DEMOCRACIES WITH FREEDOM HOUSE = 1 or 2 ##") cat("\n ############################################################## \n \n") world.file <- read.csv("FH-version.csv",header=T) alldemcy <- fullmerge[!is.na(fullmerge$fh),] alldemcy <- alldemcy[alldemcy$fh<=2,] demsbyyr <- table(alldemcy$year) } if (whichdata=="Polity") { cat("\n ##############################################################") cat("\n ## USING DEMOCRACIES WITH POLITY >= 6 ##") cat("\n ############################################################## \n \n") world.file <- read.csv("Polity-version.csv",header=T) alldemcy <- fullmerge[!is.na(fullmerge$polity),] alldemcy <- alldemcy[alldemcy$polity>=6,] demsbyyr <- table(alldemcy$year) } world <- world.file world$LNDXR <- log(world$DXR+0.3) world$ROOTGDP <- sqrt(world$GDP) world$ROOTGDPrc <- world$ROOTGDP-mean(world$ROOTGDP,na.rm=T) # re-center root GDP variable world$LNPOP <- log(world$POP) missing.rows <- (is.na(world$P) | is.na(world$LAGP) | is.na(world$ROOTGDP) | is.na(world$LNDXR) | is.na(world$LNPOP) | is.na(world$IMPORTS)) world.nomiss <- world[!missing.rows,] world.nomiss <- data.frame(country=world.nomiss$country,year=world.nomiss$year,P=world.nomiss$P,SMD=world.nomiss$SMD, LAGP=world.nomiss$LAGP,ROOTGDPrc=world.nomiss$ROOTGDPrc,LNDXR=world.nomiss$LNDXR,LNPOP=world.nomiss$LNPOP, IMPORTS=world.nomiss$IMPORTS,USAINF=world.nomiss$USAINF,OECD=world.nomiss$OECD,cnum=world.nomiss$cnum) # GDP per capita is roughly six times higher on average in the OECD than outside the OECD. # mean(world$GDP[world$year==2000 & world$OECD==1],na.rm=T) # FH: 24744; Polity 23611 # mean(world$GDP[world$year==2000 & world$OECD==0],na.rm=T) # FH: 4170; Polity 3249 ################################################################################################################# ## Number of Democracies worldwide, 1972-2000; and yearly number of SMD and PR countries yvars <- cbind(table(world$year),matrix(demsbyyr),table(world$year,world$SMD),table(world$year,world$OECD)) matplot(names(table(world$year)),yvars,xlim=c(1970,2000),type="l",lty=c(1,3,5,5,4,4),ylim=c(0,90),col=c(1,1,1,"slategray","slategray",1), cex.axis=1.3,cex.lab=1.3,las=1,lwd=2,xlab="Year",ylab="Number of Democracies") if (whichdata=="FH") { text(1989,78,"All Democracies",cex=1.3) text(1997,69,"Democracies \n in data set",cex=1.3) text(1997,47,"PR",cex=1.3) text(1997,32,"SMD",cex=1.3) text(1994,20.5,"OECD",cex=1.3) text(1990,51,"Non-OECD",cex=1.3) } #savePlot(filename=paste(whichdata,"_num_democ",sep=""),type="eps") ################################################################################################################# ## Distribution of cases among OECD and Non-OECD democracies oecd.smd<-matrix(0,nrow=3,ncol=4) oecd.smd[1:2,1:2] <- table(world$OECD,world$SMD) oecd.smd[3,1] <- sum(oecd.smd[,1]) oecd.smd[3,2] <- sum(oecd.smd[,2]) oecd.smd[1,4] <- sum(oecd.smd[1,]) oecd.smd[2,4] <- sum(oecd.smd[2,]) oecd.smd[3,4] <- sum(oecd.smd[3,]) oecd.smd[,3] <- round(100*oecd.smd[,2]/oecd.smd[,4],1) rownames(oecd.smd) <- c("Non-OECD","OECD","total") colnames(oecd.smd) <- c("PR","SMD","% SMD","total") cat("Distribution of cases among OECD and Non-OECD democracies \n \n") print(oecd.smd) ## for Freedom House data set: ## ## PR SMD % SMD total ## Non-OECD 557 458 45.1 1015 ## OECD 475 181 27.6 656 ## total 1032 639 38.2 1671 ## ## for Polity data set: ## ## PR SMD % SMD total ## Non-OECD 573 298 34.2 871 ## OECD 415 181 30.4 596 ## total 988 479 32.7 1467 ## ## Note that Polity includes 60 fewer country-years in the OECD and 144 fewer outside the OECD. ## Those 60 OECD country-years that were not rated as democratic by Polity are Spain 1977, Greece 1974, ## and all 29 years of Iceland and Luxembourg which Polity did not rate at all. ## cat("\n ************************************************************** \n \n") ################################################################################################################# ## Comparing average real prices by electoral system among OECD and non-OECD democracies, 1972-2000 oecd.pr <- matrix(0,ncol=1,nrow=29) oecd.smd <- matrix(0,ncol=1,nrow=29) nonoecd.pr <- matrix(0,ncol=1,nrow=29) nonoecd.smd <- matrix(0,ncol=1,nrow=29) for (i in 1972:2000) { oecd.pr[i-1971] <- mean(world$P[(world$year==i) & (world$OECD==1) & (world$SMD==0)],na.rm=T) oecd.smd[i-1971] <- mean(world$P[(world$year==i) & (world$OECD==1) & (world$SMD==1)],na.rm=T) nonoecd.pr[i-1971] <- mean(world$P[(world$year==i) & (world$OECD==0) & (world$SMD==0)],na.rm=T) nonoecd.smd[i-1971] <- mean(world$P[(world$year==i) & (world$OECD==0) & (world$SMD==1)],na.rm=T) } yvars <- cbind(oecd.pr,oecd.smd,nonoecd.pr,nonoecd.smd) matplot(names(table(world$year)),yvars,xlim=c(1970,2000),type="l",lty=c(1,2,4,6),col=c(1,1,1,1),ylim=c(0,140), cex.axis=1.3,cex.lab=1.3,las=1,lwd=c(2,2,2,2),xlab="Year",ylab="Real Price Index") text(1986,120,"OECD PR",cex=1.3) text(1991,95,"OECD SMD",cex=1.3) text(1979,71,"Non-OECD PR",cex=1.3) text(1979,42,"Non-OECD SMD",cex=1.3) #savePlot(filename=paste(whichdata,"_avg_prices",sep=""),type="eps") ################################################################################################################# ## Data missingness in world sample. miss <- data.frame(VARIABLES=c("Real prices","Imports as percent of GDP","Population","GDP per capita", "Yearly change in exchange rate","United States inflation rate","SMD electoral rules","Overall"), MISSINGNESS=round(c(((table(is.na(world$P))[2])/nrow(world)),((table(is.na(world$IMPORTS))[2])/nrow(world)), ((table(is.na(world$POP))[2])/nrow(world)),((table(is.na(world$GDP))[2])/nrow(world)), ((table(is.na(world$LNDXR))[2])/nrow(world)),((table(is.na(world$USAINF))[2])/nrow(world)), ((table(is.na(world$SMD))[2])/nrow(world)),(sum(is.na(cbind(world$P,world$IMPORTS,world$POP, world$GDP,world$LNDXR,world$USAINF,world$SMD)))/(7*nrow(world)))),5)) cat("Data missingness \n \n") print(miss) ## for Freedom House data set: ## ## VARIABLES MISSINGNESS ## 1 Real prices 0.13166 ## 2 Imports as percent of GDP 0.06284 ## 3 Population 0.04309 ## 4 GDP per capita 0.04010 ## 5 Yearly change in exchange rate 0.03052 ## 6 United States inflation rate NA ## 7 SMD electoral rules NA ## 8 Overall 0.04403 ## ## For Polity data set: ## ## VARIABLES MISSINGNESS ## 1 Real prices 0.03067 ## 2 Imports as percent of GDP 0.01977 ## 3 Population NA ## 4 GDP per capita 0.00273 ## 5 Yearly change in exchange rate 0.01704 ## 6 United States inflation rate NA ## 7 SMD electoral rules NA ## 8 Overall 0.01003 ## cat("\n ************************************************************** \n \n") ################################################################################################################# ## Scatterplot of Real Prices vs. per capita GDP, with linear fits ## for square-root transformed and log-transformed per capita GDP. p.oecd <- world$P p.oecd[world$OECD==0] <- NA p.non <- world$P p.non[world$OECD==1] <- NA yvars <- cbind(p.oecd,p.non) matplot(world$GDP,yvars,type="p",pch=c(1,3),col=c(1,1),xlim=c(0,50000),ylim=c(0,200), cex.axis=1.3,cex.lab=1.3,las=1,lwd=c(1,1),xlab="Per capita GDP",ylab="Real Price Index") rootfit <- lm(world$P~sqrt(world$GDP)) logfit <- lm(world$P~log(world$GDP)) rootline <- matrix(NA,nrow=460,ncol=2) logline <- matrix(NA,nrow=460,ncol=2) for (i in 1:460) { rootline[i,] <- c(i*100,coefficients(rootfit)[1] + coefficients(rootfit)[2]*sqrt(i*100)) logline[i,] <- c(i*100,coefficients(logfit)[1] + coefficients(logfit)[2]*log(i*100)) } points(rootline,type="l",lwd=3,lty=1); text(45500,164,"Square Root",cex=1.2) points(logline,type="l",lwd=3,lty=2); text(45500,131,"Log",cex=1.2) #savePlot(filename=paste(whichdata,"_prices_GDP",sep=""),type="eps") ################################################################################################################# ## Some analyses of the relationship between real price levels (P) and ## square-root transformed GDP per capita (ROOTGDP), as pooled vs. panel data, to ## compare and illustrate the similarities/differences between OLS and linear GEE. plot(world$ROOTGDP,world$P,cex=0.5,xlim=c(0,220),ylim=c(0,200),cex.axis=1.3,cex.lab=1.3, las=1,xlab="Per capita GDP, square root",ylab="Real Price Index",col="darkgray") for (i in 1:max(world$cnum)) { useobs <- (world$cnum == i) if ((sum(!is.na(world$P[useobs])) > 3) & (sum(!is.na(world$ROOTGDP[useobs])) > 3)) { ols <- lm(P~ROOTGDP,world[useobs,]) #cat(i,"\n") #abline(coefficients(ols),lty=3) #plot(world$ROOTGDP,world$P,cex=1,xlim=c(0,220),ylim=c(0,200),cex.axis=1.3,cex.lab=1.3,las=1,xlab="Per capita GDP, square root",ylab="Real Price Index",pch=".") #points(world$ROOTGDP[useobs],world$P[useobs],cex=2,lwd=2) minx <- min(world$ROOTGDP[useobs][!is.na(world$P[useobs])],na.rm=T) maxx <- max(world$ROOTGDP[useobs][!is.na(world$P[useobs])],na.rm=T) lines(c(minx,maxx),c(coefficients(ols)[1]+(minx*coefficients(ols)[2]),coefficients(ols)[1]+(maxx*coefficients(ols)[2])),lwd=0.5) } } minx <- min(world$ROOTGDP[!is.na(world$P)],na.rm=T) maxx <- max(world$ROOTGDP[!is.na(world$P)],na.rm=T) gee.ind.bivar <- gee(P~ROOTGDP,id=cnum,data=world,corstr="independence") # equivalent to pooled OLS lines(c(minx,maxx),c(coefficients(gee.ind.bivar)[1]+(minx*coefficients(gee.ind.bivar)[2]), coefficients(gee.ind.bivar)[1]+(maxx*coefficients(gee.ind.bivar)[2])),lwd=3,lty=4) gee.exch.bivar <- gee(P~ROOTGDP,id=cnum,data=world,corstr="exchangeable") lines(c(minx,maxx),c(coefficients(gee.exch.bivar)[1]+(minx*coefficients(gee.exch.bivar)[2]), coefficients(gee.exch.bivar)[1]+(maxx*coefficients(gee.exch.bivar)[2])),lwd=3,lty=1) text(215,145,"Pooled\n OLS",cex=1.1) text(180,85,"GEE\n Exchangeable",cex=1.1) #savePlot(filename=paste(whichdata,"_prices_GDP_root_fits",sep=""),type="eps") ################################################################################################################# ## Show that the four control variables are distributed similarly between SMD and PR countries. ## SMD density plots are denoted with solid lines; PR density plots denoted with dashed lines. par(mar=c(4,4,2,1),mfrow=c(2,2)) world.SMD <- world[world$SMD==1,] world.PR <- world[world$SMD==0,] plot(density(world.SMD$ROOTGDP,na.rm=T),lwd=2,lty=1,main="per capita GDP, square root",xlab="") lines(density(world.PR$ROOTGDP,na.rm=T),lwd=2,lty=2) plot(density(world.SMD$LNDXR,na.rm=T),lwd=2,lty=1,main="Yearly change in exchange rates, log",xlab="") lines(density(world.PR$LNDXR,na.rm=T),lwd=2,lty=2) plot(density(world.PR$LNPOP,na.rm=T),lwd=2,lty=2,main="Population, log",xlab="",xlim=c(8,22)) lines(density(world.SMD$LNPOP,na.rm=T),lwd=2,lty=1) plot(density(world.PR$IMPORTS,na.rm=T),lwd=2,lty=2,main="Imports as a percentage of GDP",xlab="") lines(density(world.SMD$IMPORTS,na.rm=T),lwd=2,lty=1) #savePlot(filename=paste(whichdata,"_matched_IVs",sep=""),type="eps") qqplot(world.SMD$ROOTGDP,world.PR$ROOTGDP,main="per capita GDP, square root",xlab="majoritarian",ylab="proportional") abline(0,1) qqplot(world.SMD$LNDXR,world.PR$LNDXR,main="Yearly change in exchange rates, log",xlab="majoritarian",ylab="proportional") # worst one abline(0,1) qqplot(world.SMD$LNPOP,world.PR$LNPOP,main="Population, log",xlab="majoritarian",ylab="proportional") abline(0,1) qqplot(world.SMD$IMPORTS,world.PR$IMPORTS,main="Imports as a percentage of GDP",xlab="majoritarian",ylab="proportional") abline(0,1) #savePlot(filename=paste(whichdata,"_matched_IVs_qqplot",sep=""),type="eps") par(mar=c(5,4.5,0.5,0.5),mfrow=c(1,1)) ################################################################################################################# ## GEE Models. Dependent variable is real prices. ## For documentation, see http://finzi.psych.upenn.edu/R/library/gee/html/gee.html p.listwise.del <- world$P[!missing.rows] ## ## MODEL: GEE Exchangeable ## cat("\n ******************************************** ") cat("\n ********* GEE Model: Exchangeable ********** \n \n") GEE.exch <- gee(P~SMD+LAGP+ROOTGDP+LNDXR+LNPOP+IMPORTS+USAINF, id=cnum, data=world, corstr="exchangeable") print(summary(GEE.exch)$coefficients) long.term <- 100*(GEE.exch$coefficients[2]/(1-GEE.exch$coefficients[3]))/mean(world$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") plot(p.listwise.del,GEE.exch$fitted.values,cex.axis=1.3,cex.lab=1.3,xlab="Observed real price levels",ylab="Predicted real price levels: GEE exchangeable") abline(0,1) n <- c(names(GEE.exch$coefficients),"SMDxROOTGDP") outtable <- data.frame(label=c(matrix(sapply(n,c," SE")),"","Observations","rho","","Long-term SMD effect")) outtable$GEE.exch <- c(round(matrix(t(summary(GEE.exch)$coefficients[,c(1,4)])),3), # coefficients and SEs "","","",GEE.exch$nobs, # number of observations round(GEE.exch$working.correlation[1,2],3),"", # correlation parameter (rho) round(long.term,1)) # long term effect of SMD ## ## MODEL: GEE Exchangeable with GDP interaction term ## cat("\n ******************************************************* ") cat("\n ****** GEE Model: Exchangeable + GDP Interaction ****** \n \n") GEE.exch.inter <- gee(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc, id=cnum, data=world, corstr="exchangeable") print(summary(GEE.exch.inter)$coefficients) cat("\n Variance of residuals:",var(GEE.exch$residuals),"\n \n") long.term <- 100*(GEE.exch.inter$coefficients[2]/(1-GEE.exch.inter$coefficients[3]))/mean(world$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") outtable$GEE.exch.inter <- c(round(matrix(t(summary(GEE.exch.inter)$coefficients[,c(1,4)])),3), # coefficients and SEs "",GEE.exch.inter$nobs, # number of observations round(GEE.exch.inter$working.correlation[1,2],3),"", # correlation parameter (rho) round(long.term,1)) # long term effect of SMD at mean ROOTGDP yhat <- data.frame(GDPpc=rep(NA,204),SMD=NA,PR=NA,SMD_PR=NA) pos <- 1 for (i in -65:138) { yhat$GDPpc[pos] <- (i+mean(world$ROOTGDP,na.rm=T))^2 iv.smd <- c(1,1,mean(world$LAGP,na.rm=T),i,mean(world$LNDXR,na.rm=T),mean(world$LNPOP,na.rm=T),mean(world$IMPORTS,na.rm=T),mean(world$USAINF,na.rm=T),i) iv.pr <- c(1,0,mean(world$LAGP,na.rm=T),i,mean(world$LNDXR,na.rm=T),mean(world$LNPOP,na.rm=T),mean(world$IMPORTS,na.rm=T),mean(world$USAINF,na.rm=T),0) yhat$SMD[pos] <- iv.smd %*% GEE.exch.inter$coefficients yhat$PR[pos] <- iv.pr %*% GEE.exch.inter$coefficients pos <- pos+1 } yhat$SMD_PR <- yhat$SMD-yhat$PR matplot(yhat$GDPpc,cbind(yhat$SMD,yhat$PR),ylim=c(55,95),type="l",lwd=2,col=1,cex.axis=1.3,cex.lab=1.3,xlab="per capita GDP",ylab="Predicted real price level") #savePlot(filename=paste(whichdata,"_pv_alldems",sep=""),type="eps") ## ## MODEL: GEE Exchangeable with OECD interaction term ## cat("\n ************************************************************** ") cat("\n ****** GEE Model: Exchangeable + OECD Interaction ****** \n \n") GEE.exch.oecd <- gee(P~SMD+LAGP+OECD+LNDXR+LNPOP+IMPORTS+USAINF+SMD:OECD, id=cnum, data=world, corstr="exchangeable") print(summary(GEE.exch.oecd)$coefficients) long.term <- 100*(GEE.exch.oecd$coefficients[2]/(1-GEE.exch.oecd$coefficients[3]))/mean(world$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") outtable$GEE.exch.oecd <- c(round(matrix(t(summary(GEE.exch.oecd)$coefficients[,c(1,4)])),3), # coefficients and SEs "",GEE.exch.oecd$nobs, # number of observations round(GEE.exch.oecd$working.correlation[1,2],3),"", # correlation parameter (rho) round(long.term,1)) # long term effect of SMD at mean ROOTGDP iv.mat <- cbind(c(1,1,mean(world$LAGP,na.rm=T),1,mean(world$LNDXR,na.rm=T),mean(world$LNPOP,na.rm=T),mean(world$IMPORTS,na.rm=T),mean(world$USAINF,na.rm=T),1), c(1,0,mean(world$LAGP,na.rm=T),1,mean(world$LNDXR,na.rm=T),mean(world$LNPOP,na.rm=T),mean(world$IMPORTS,na.rm=T),mean(world$USAINF,na.rm=T),0), c(1,1,mean(world$LAGP,na.rm=T),0,mean(world$LNDXR,na.rm=T),mean(world$LNPOP,na.rm=T),mean(world$IMPORTS,na.rm=T),mean(world$USAINF,na.rm=T),0), c(1,0,mean(world$LAGP,na.rm=T),0,mean(world$LNDXR,na.rm=T),mean(world$LNPOP,na.rm=T),mean(world$IMPORTS,na.rm=T),mean(world$USAINF,na.rm=T),0)) yhat <- GEE.exch.oecd$coefficients %*% iv.mat colnames(yhat) <- c("OECD-SMD","OECD-PR","non-SMD","non-PR") print(t(yhat)) cat("OECD: SMD-PR =",yhat[1]-yhat[2],"\n") cat("non-OECD: SMD-PR =",yhat[3]-yhat[4],"\n") ## ## MODEL: GEE Exchangeable with GDP interaction term and control for US account balance rather than inflation rate ## cat("\n ******************************************************************************** ") cat("\n ****** GEE Model: Exchangeable + Interaction + US account balance control ****** \n \n") GEE.exch.usab <- gee(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAB+SMD:ROOTGDPrc, id=cnum, data=world, corstr="exchangeable") print(summary(GEE.exch.usab)$coefficients) long.term <- 100*(GEE.exch.usab$coefficients[2]/(1-GEE.exch.usab$coefficients[3]))/mean(world$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") outtable$GEE.exch.usab <- c(round(matrix(t(summary(GEE.exch.usab)$coefficients[,c(1,4)])),3), # coefficients and SEs "",GEE.exch.usab$nobs, # number of observations round(GEE.exch.usab$working.correlation[1,2],3),"", # correlation parameter (rho) round(long.term,1)) # long term effect of SMD at mean ROOTGDP ## ## MODEL: GEE Exchangeable excluding countries with mixed electoral systems ## This excludes 115 country-years in the Freedom House data set and 152 years in the Polity data set. ## cat("\n **************************************************************** ") cat("\n ********* GEE Model: Exchangeable, excluding mixed ES ********** \n \n") world.nomix <- world[world$MIX==0,] world.nomix$ROOTGDPrc <- world.nomix$ROOTGDP-mean(world.nomix$ROOTGDP,na.rm=T) # re-center root GDP variable GEE.exch.nomix <- gee(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc, id=cnum, data=world.nomix, corstr="exchangeable") print(summary(GEE.exch.nomix)$coefficients) long.term <- 100*(GEE.exch.nomix$coefficients[2]/(1-GEE.exch.nomix$coefficients[3]))/mean(world.nomix$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") outtable$GEE.exch.nomix <- c(round(matrix(t(summary(GEE.exch.nomix)$coefficients[,c(1,4)])),3), # coefficients and SEs "",GEE.exch.nomix$nobs, # number of observations round(GEE.exch.nomix$working.correlation[1,2],3),"", # correlation parameter (rho) round(long.term,1)) # long term effect of SMD yhat <- data.frame(GDPpc=rep(NA,204),SMD=NA,PR=NA,SMD_PR=NA) pos <- 1 for (i in -65:138) { yhat$GDPpc[pos] <- (i+mean(world.nomix$ROOTGDP,na.rm=T))^2 iv.smd <- c(1,1,mean(world.nomix$LAGP,na.rm=T),i,mean(world.nomix$LNDXR,na.rm=T),mean(world.nomix$LNPOP,na.rm=T),mean(world.nomix$IMPORTS,na.rm=T),mean(world.nomix$USAINF,na.rm=T),i) iv.pr <- c(1,0,mean(world.nomix$LAGP,na.rm=T),i,mean(world.nomix$LNDXR,na.rm=T),mean(world.nomix$LNPOP,na.rm=T),mean(world.nomix$IMPORTS,na.rm=T),mean(world.nomix$USAINF,na.rm=T),0) yhat$SMD[pos] <- iv.smd %*% GEE.exch.nomix$coefficients yhat$PR[pos] <- iv.pr %*% GEE.exch.nomix$coefficients pos <- pos+1 } yhat$SMD_PR <- yhat$SMD-yhat$PR matplot(yhat$GDPpc,cbind(yhat$SMD,yhat$PR),ylim=c(55,95),type="l",lwd=2,col=1,cex.axis=1.3,cex.lab=1.3,xlab="per capita GDP (no mixed ES)",ylab="Predicted real price level") #savePlot(filename=paste(whichdata,"_pv_nomixedES",sep=""),type="eps") ## ## MODEL: GEE Exchangeable with all Japan and Mongolia observations coded as SMD=1 ## cat("\n ******************************************************************************** ") cat("\n ********* GEE Model: Exchangeable, recoding Japan and Mongolia as SMD ********** \n \n") world.jpnmng <- world world.jpnmng$SMD[world$WBcode=="JPN"] <- 1 world.jpnmng$SMD[world$WBcode=="MNG"] <- 1 GEE.exch.jpnmng <- gee(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc, id=cnum, data=world.jpnmng, corstr="exchangeable") print(summary(GEE.exch.jpnmng)$coefficients) long.term <- 100*(GEE.exch.jpnmng$coefficients[2]/(1-GEE.exch.jpnmng$coefficients[3]))/mean(world$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") outtable$GEE.exch.jpnmng <- c(round(matrix(t(summary(GEE.exch.jpnmng)$coefficients[,c(1,4)])),3), # coefficients and SEs "",GEE.exch.jpnmng$nobs, # number of observations round(GEE.exch.jpnmng$working.correlation[1,2],3),"", # correlation parameter (rho) round(long.term,1)) # long term effect of SMD ## ## MODEL: GEE AR(1) ## cat("\n ******************************************** ") cat("\n ************ GEE Model: AR(1) ************** \n \n") GEE.AR1 <- gee(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc, id=cnum, data=world, corstr="AR-M", Mv=1) print(summary(GEE.AR1)$coefficients) cat("\n Variance of residuals:",var(GEE.AR1$residuals),"\n \n") long.term <- 100*(GEE.AR1$coefficients[2]/(1-GEE.AR1$coefficients[3]))/mean(world$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") plot(p.listwise.del,GEE.AR1$fitted.values,cex.axis=1.3,cex.lab=1.3,xlab="Observed real price levels",ylab="Predicted real price levels: GEE AR(1)") abline(0,1) outtable$GEE.AR1 <- c(round(matrix(t(summary(GEE.AR1)$coefficients[,c(1,4)])),3), # coefficients and SEs "",GEE.AR1$nobs, # number of observations round(GEE.AR1$working.correlation[1,2],3),"", # correlation parameter (rho) round(long.term,1)) # long term effect of SMD ## ## MODEL: GEE Exchangeability, imputed datasets ## cat("\n *************************************************** ") cat("\n ******** GEE Model: Exchangeable + Imputed ******** \n \n") ################################################################################################################# ## Imputation using Amelia II preimpute <- world preimpute$LNLANDPC <- log(preimpute$LANDPC) preimpute$cy <- preimpute$country <- preimpute$XR <- preimpute$DXR <- NULL preimpute$PPP <- preimpute$POP <- preimpute$LANDPC <- preimpute$GDP <- NULL preimpute$LAGP <- preimpute$ROOTGDPrc <- NULL # note: include variables such as ETHNIC and HONGOV in the imputation model for better guesses. world.imp <- amelia(preimpute,m=10,idvars="WBcode",ts="year",cs="cnum",lags="P",polytime=1,noms=c("SMD","OECD","MIX"),write.out=F) # to save copies of the imputed data sets to disk, change the 'write.out' option to TRUE. # note: the 'intercs' option does not improve the imputation, is computationally slow with # 'polytime'>1, and requires a large 'empri' (~160) as well as setting 'incheck'=F. createlagP <- function(imp,orig) { ## this function creates a lagged P variable that is equal to the actual previous year's P ## if observed, or the imputed previous year's P if missing. In a small number of cases, ## the previous year's P was neither observed nor imputed because that year was not in the ## data set. In those cases, the lagged P variable remains coded as missing (NA). ret <- imp ret$LAGP <- orig$LAGP for (i in 2:nrow(ret)) { if ((is.na(ret$LAGP[i])) & (ret$year[i]==(ret$year[i-1]+1)) & (ret$WBcode[i]==ret$WBcode[i-1])) { ret$LAGP[i] <- ret$P[i-1] } } ret } for (i in 1:10) { world.imp[[i]] <- createlagP(world.imp[[i]],world) # replace missing lag P values with imputed values world.imp[[i]] <- transform(world.imp[[i]],ROOTGDPrc=(ROOTGDP-mean(ROOTGDP,na.rm=T))) } ################################################################################################################# ## Results of multiple imputation: mean values of five imputed variables across 10 imputed datasets mean.vals <- matrix(NA,nrow=11,ncol=5) mean.vals[1,] <- rbind(mean(world$P,na.rm=T),mean(world$IMPORTS,na.rm=T),mean(world$LNPOP,na.rm=T),mean(world$ROOTGDP,na.rm=T),mean(world$LNDXR,na.rm=T)) for (i in 1:10) { mean.vals[(i+1),] <- rbind(mean(world.imp[[i]]$P),mean(world.imp[[i]]$IMPORTS),mean(world.imp[[i]]$LNPOP), mean(world.imp[[i]]$ROOTGDP),mean(world.imp[[i]]$LNDXR)) } cat("Mean values of five imputed variables across 10 imputed datasets \n \n") colnames(mean.vals) <- c("P","IMPORTS","LNPOP","ROOTGDP","LNDXR") rownames(mean.vals) <- c("Original sample","Imputed 1","Imputed 2","Imputed 3","Imputed 4","Imputed 5","Imputed 6","Imputed 7","Imputed 8","Imputed 9","Imputed 10") print(mean.vals) cat("\n ************************************************************** \n \n") ################################################################################################################# ## Fit GEE models to imputed datasets and perform model averaging GEE.imp <- list() long.term <- NULL for(i in 1:10) { cat("Estimate GEE model, imputed dataset",i,"\n") GEE.imp[[i]] <- gee(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc, id=cnum, data=world.imp[[i]], corstr="exchangeable") long.term <- c(long.term,100*(GEE.imp[[i]]$coefficients[2]/(1-GEE.imp[[i]]$coefficients[3]))/mean(world.imp[[i]]$P,na.rm=T)) } coeff.imp <- cbind(coefficients(GEE.imp[[1]]),coefficients(GEE.imp[[2]]),coefficients(GEE.imp[[3]]),coefficients(GEE.imp[[4]]),coefficients(GEE.imp[[5]]), coefficients(GEE.imp[[6]]),coefficients(GEE.imp[[7]]),coefficients(GEE.imp[[8]]),coefficients(GEE.imp[[9]]),coefficients(GEE.imp[[10]])) coeff.mean <- rowMeans(coeff.imp) naivese.imp <- cbind(summary(GEE.imp[[1]])$coefficients[,2],summary(GEE.imp[[2]])$coefficients[,2],summary(GEE.imp[[3]])$coefficients[,2], summary(GEE.imp[[4]])$coefficients[,2],summary(GEE.imp[[5]])$coefficients[,2],summary(GEE.imp[[6]])$coefficients[,2], summary(GEE.imp[[7]])$coefficients[,2],summary(GEE.imp[[8]])$coefficients[,2],summary(GEE.imp[[9]])$coefficients[,2], summary(GEE.imp[[10]])$coefficients[,2]) robustse.imp <- cbind(summary(GEE.imp[[1]])$coefficients[,4],summary(GEE.imp[[2]])$coefficients[,4],summary(GEE.imp[[3]])$coefficients[,4], summary(GEE.imp[[4]])$coefficients[,4],summary(GEE.imp[[5]])$coefficients[,4],summary(GEE.imp[[6]])$coefficients[,4], summary(GEE.imp[[7]])$coefficients[,4],summary(GEE.imp[[8]])$coefficients[,4],summary(GEE.imp[[9]])$coefficients[,4], summary(GEE.imp[[10]])$coefficients[,4]) Sq.squared <- 10*rowMeans((coeff.imp-coeff.mean)^2)/9 naivese.mean <- sqrt(rowMeans(naivese.imp^2) + (Sq.squared*(1+(1/10)))) robustse.mean <- sqrt(rowMeans(robustse.imp^2) + (Sq.squared*(1+(1/10)))) correl.mean <- mean(c(GEE.imp[[1]]$working.correlation[1,2],GEE.imp[[2]]$working.correlation[1,2],GEE.imp[[3]]$working.correlation[1,2], GEE.imp[[4]]$working.correlation[1,2],GEE.imp[[5]]$working.correlation[1,2],GEE.imp[[6]]$working.correlation[1,2], GEE.imp[[7]]$working.correlation[1,2],GEE.imp[[8]]$working.correlation[1,2],GEE.imp[[9]]$working.correlation[1,2], GEE.imp[[10]]$working.correlation[1,2])) long.term <- mean(long.term) cat("\n \n Results from model averaging across 10 imputed datasets: \n \n") print(rbind(coeff.mean,naivese.mean,robustse.mean)) cat("\n") cat("rho=",correl.mean,"\n") cat("N=",GEE.imp[[1]]$nobs,"\n") cat("Long-term effect of switch to SMD:",long.term,"% \n \n") outtable$GEE.imputed <- c(round(matrix(t(cbind(coeff.mean,robustse.mean))),3), # coefficients and SEs "",GEE.imp[[1]]$nobs, # number of observations round(correl.mean,3),"", # correlation parameter (rho) round(long.term,1)) # long term effect of SMD ## ## MODEL: OLS with PCSE, with and without country fixed effects ## cat("\n **************************************************** ") cat("\n ******************** Model: OLS ******************** \n \n") OLS <- lm(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc,data=world.nomiss) print(summary(OLS)$coefficients) cat("Note: these are regular standard errors, not PCSEs \n") long.term <- 100*(OLS$coefficients[2]/(1-OLS$coefficients[3]))/mean(world$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") outtable$OLS <- c(round(matrix(t(cbind(summary(OLS)$coefficients[,c(1,2)]))),3), # coefficients and SEs "",length(OLS$residuals), # number of observations round(summary(OLS)$r.squared,3),"", # R^2 of model round(long.term,1)) # long term effect of SMD if (whichdata=="FH") { cfe <- data.frame(C1=rep(0,nrow(world)),C2=0,C3=0,C4=0,C5=0,C6=0,C7=0,C8=0,C9=0,C10=0,C11=0, C12=0,C13=0,C14=0,C15=0,C16=0,C17=0,C18=0,C19=0,C20=0,C21=0,C22=0,C23=0, C24=0,C25=0,C26=0,C27=0,C28=0,C29=0,C30=0,C31=0,C32=0,C33=0,C34=0,C35=0, C36=0,C37=0,C38=0,C39=0,C40=0,C41=0,C42=0,C43=0,C44=0,C45=0,C46=0,C47=0, C48=0,C49=0,C50=0,C51=0,C52=0,C53=0,C54=0,C55=0,C56=0,C57=0,C58=0,C59=0, C60=0,C61=0,C62=0,C63=0,C64=0,C65=0,C66=0,C67=0,C68=0,C69=0,C70=0,C71=0, C72=0,C73=0,C74=0,C75=0,C76=0,C77=0,C78=0,C79=0,C80=0,C81=0,C82=0,C83=0, C84=0,C85=0,C86=0,C87=0,C88=0,C89=0,C90=0) } else { cfe <- data.frame(C1=rep(0,nrow(world)),C2=0,C3=0,C4=0,C5=0,C6=0,C7=0,C8=0,C9=0,C10=0,C11=0, C12=0,C13=0,C14=0,C15=0,C16=0,C17=0,C18=0,C19=0,C20=0,C21=0,C22=0,C23=0, C24=0,C25=0,C26=0,C27=0,C28=0,C29=0,C30=0,C31=0,C32=0,C33=0,C34=0,C35=0, C36=0,C37=0,C38=0,C39=0,C40=0,C41=0,C42=0,C43=0,C44=0,C45=0,C46=0,C47=0, C48=0,C49=0,C50=0,C51=0,C52=0,C53=0,C54=0,C55=0,C56=0,C57=0,C58=0,C59=0, C60=0,C61=0,C62=0,C63=0,C64=0,C65=0,C66=0,C67=0,C68=0,C69=0,C70=0,C71=0, C72=0,C73=0,C74=0,C75=0,C76=0) } for (i in 1:max(world$cnum)) cfe[world$cnum==i,i] <- 1 world.cfe <- world world.cfe$GDPSMD <- world.cfe$SMD * world.cfe$ROOTGDPrc world.cfe <- cbind(world.cfe,cfe) #write.csv(world.cfe,paste(whichdata,"_cfe.csv",sep=""),na="",row.names=F) ## ## Need to get PCSEs in STATA. (These results are from STATA 7.) ## ## Open the (FH/Polity)_cfe.csv data file in STATA and enter the commands: ## ## tsset cnum year ## set matsize 250 ## xtpcse p smd lagp rootgdprc lndxr lnpop imports usainf gdpsmd, pa ## ## For the Freedom House data set, the results are: ## ## Number of gaps in sample: 14 ## (note: at least one disturbance covariance assumed 0, no common time periods ## between panels) ## ## Linear regression, correlated panels corrected standard errors (PCSEs) ## ## Group variable: cnum Number of obs = 1410 ## Time variable: year Number of groups = 76 ## Panels: correlated (unbalanced) Obs per group: min = 4 ## Autocorrelation: no autocorrelation avg = 18.55263 ## Sigma computed by pairwise selection max = 29 ## Estimated covariances = 2926 R-squared = 0.9589 ## Estimated autocorrelations = 0 Wald chi2(8) = 9983.90 ## Estimated coefficients = 9 Prob > chi2 = 0.0000 ## ## ------------------------------------------------------------------------------ ## | Panel-corrected ## | Coef. Std. Err. z P>|z| [95% Conf. Interval] ## -------------+---------------------------------------------------------------- ## smd | -2.772249 .5359559 -5.17 0.000 -3.822704 -1.721795 ## lagp | .8285367 .0346388 23.92 0.000 .7606459 .8964274 ## rootgdprc | .0880011 .0304248 2.89 0.004 .0283695 .1476327 ## lndxr | -7.321363 .9932743 -7.37 0.000 -9.268144 -5.374581 ## lnpop | -.5506087 .1276209 -4.31 0.000 -.8007412 -.3004763 ## imports | -.1038228 .0180968 -5.74 0.000 -.1392919 -.0683537 ## usainf | .5646281 .1871105 3.02 0.003 .1978982 .9313579 ## gdpsmd | -.0074393 .0129549 -0.57 0.566 -.0328304 .0179519 ## _cons | 16.34105 4.095762 3.99 0.000 8.313507 24.3686 ## ------------------------------------------------------------------------------ ## ## ## For the Polity IV data set, the results are: ## ## Number of gaps in sample: 5 ## (note: at least one disturbance covariance assumed 0, no common time periods ## between panels) ## ## Linear regression, correlated panels corrected standard errors (PCSEs) ## ## Group variable: cnum Number of obs = 1389 ## Time variable: year Number of groups = 74 ## Panels: correlated (unbalanced) Obs per group: min = 5 ## Autocorrelation: no autocorrelation avg = 18.77027 ## Sigma computed by pairwise selection max = 29 ## Estimated covariances = 2775 R-squared = 0.9613 ## Estimated autocorrelations = 0 Wald chi2(8) = 9221.25 ## Estimated coefficients = 9 Prob > chi2 = 0.0000 ## ## ------------------------------------------------------------------------------ ## | Panel-corrected ## | Coef. Std. Err. z P>|z| [95% Conf. Interval] ## -------------+---------------------------------------------------------------- ## smd | -2.389164 .523514 -4.56 0.000 -3.415233 -1.363096 ## lagp | .7921301 .0371786 21.31 0.000 .7192615 .8649988 ## rootgdprc | .1226993 .0320867 3.82 0.000 .0598105 .1855881 ## lndxr | -5.795914 .8665176 -6.69 0.000 -7.494258 -4.097571 ## lnpop | -.5106625 .152031 -3.36 0.001 -.8086378 -.2126872 ## imports | -.1049667 .0199693 -5.26 0.000 -.1441058 -.0658275 ## usainf | .7105872 .2071833 3.43 0.001 .3045155 1.116659 ## gdpsmd | -.0154128 .0113002 -1.36 0.173 -.0375608 .0067352 ## _cons | 18.51463 4.808991 3.85 0.000 9.089177 27.94008 ## ------------------------------------------------------------------------------ ## ## ## Note that it is pointless to fit this model with country fixed effects because ## of the nature of the IVs SMD and LNPOP. SMD is nearly collinear with the FEs, ## and LNPOP is extremely sluggish over time. ## cat("\n ***************************************************** ") cat("\n ************ Model: OLS with country FEs ************ \n \n") if (whichdata=="FH") { OLS.cfe <- lm(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc+C1+C2+C3+ C4+C5+C6+C7+C8+C9+C10+C11+C12+C13+C14+C15+C16+C17+C18+C19+C20+C21+C22+ C23+C24+C25+C26+C27+C28+C29+C30+C31+C32+C33+C34+C35+C36+C37+C38+C39+ C40+C41+C42+C43+C44+C45+C46+C47+C48+C49+C50+C51+C52+C53+C54+C55+C56+ C57+C58+C59+C60+C61+C62+C63+C64+C65+C66+C67+C68+C69+C70+C71+C72+C73+ C74+C75+C76+C77+C78+C79+C80+C81+C82+C83+C84+C85+C86+C88+C89+C90, data=world.cfe,na.action=na.omit) } else { OLS.cfe <- lm(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc+C1+C2+C3+ C4+C5+C6+C7+C8+C9+C10+C11+C12+C13+C14+C15+C16+C17+C18+C19+C20+C21+C22+ C23+C24+C25+C26+C27+C28+C29+C30+C31+C32+C33+C34+C35+C36+C37+C38+C39+ C40+C41+C42+C43+C44+C45+C46+C47+C48+C49+C50+C51+C52+C53+C54+C55+C56+ C57+C58+C59+C60+C61+C62+C63+C64+C65+C66+C67+C68+C69+C70+C71+C72+C73+ C74+C75+C76,data=world.cfe,na.action=na.omit) } print(summary(OLS.cfe)$coefficients[c(1:8,nrow(summary(OLS.cfe)$coefficients)),]) cat("Note: these are regular standard errors, not PCSEs \n") cat("fixed effects not shown \n") ## ## To get the PCSEs in STATA, for the FH data set, use the command: ## xtpcse p smd lagp rootgdprc lndxr lnpop imports usainf gdpsmd c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20 c21 c22 c23 c24 c25 c26 c27 c28 c29 c30 c31 c32 c33 c34 c35 c36 c37 c38 c39 c40 c41 c42 c43 c44 c45 c46 c47 c48 c49 c50 c51 c52 c53 c54 c55 c56 c57 c58 c59 c60 c61 c62 c63 c64 c65 c66 c67 c68 c69 c70 c71 c72 c73 c74 c75 c76 c77 c78 c79 c80 c81 c82 c83 c84 c85 c86 c88 c89 c90, pa ## ## ------------------------------------------------------------------------------ ## | Panel-corrected ## | Coef. Std. Err. z P>|z| [95% Conf. Interval] ## -------------+---------------------------------------------------------------- ## smd | 1.897621 1.050417 1.81 0.071 -.1611596 3.956401 ## lagp | .69138 .046957 14.72 0.000 .599346 .783414 ## rootgdprc | .1149294 .0382921 3.00 0.003 .0398782 .1899806 ## lndxr | -9.041725 .9744763 -9.28 0.000 -10.95166 -7.131786 ## lnpop | -4.963701 3.304244 -1.50 0.133 -11.4399 1.512497 ## imports | -.1981682 .0342005 -5.79 0.000 -.2652 -.1311364 ## usainf | .6750369 .1875908 3.60 0.000 .3073657 1.042708 ## gdpsmd | -.0135296 .0160291 -0.84 0.399 -.044946 .0178867 ## _cons | 99.3294 57.6406 1.72 0.085 -13.6441 212.3029 ## ## And for the Polity IV data set, use the command: ## xtpcse p smd lagp rootgdprc lndxr lnpop imports usainf gdpsmd c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20 c21 c22 c23 c24 c25 c26 c27 c28 c29 c30 c31 c32 c33 c34 c35 c36 c37 c38 c39 c40 c41 c42 c43 c44 c45 c46 c47 c48 c49 c50 c51 c52 c53 c54 c55 c56 c57 c58 c59 c60 c61 c62 c63 c64 c65 c66 c67 c68 c69 c70 c71 c72 c73 c74 c75 c76, pa ## ## ------------------------------------------------------------------------------ ## | Panel-corrected ## | Coef. Std. Err. z P>|z| [95% Conf. Interval] ## -------------+---------------------------------------------------------------- ## smd | 1.454401 .9602112 1.51 0.130 -.4275786 3.33638 ## lagp | .6377219 .0466521 13.67 0.000 .5462855 .7291584 ## rootgdprc | .160353 .0397638 4.03 0.000 .0824174 .2382885 ## lndxr | -6.79972 .8012776 -8.49 0.000 -8.370195 -5.229245 ## lnpop | -5.418723 2.944589 -1.84 0.066 -11.19001 .3525657 ## imports | -.248796 .0375677 -6.62 0.000 -.3224275 -.1751646 ## usainf | .8265246 .2098783 3.94 0.000 .4151707 1.237879 ## gdpsmd | -.0279008 .0159547 -1.75 0.080 -.0591715 .0033698 ## _cons | 110.6734 42.30833 2.62 0.009 27.75055 193.5962 ## ## The coefficient on LNPOP has ballooned from less than one to -5. ## Hence the fact that the coefficient on SMD loses its significance is substantively uninformative. ## long.term <- 100*(OLS.cfe$coefficients[2]/(1-OLS.cfe$coefficients[3]))/mean(world.cfe$P,na.rm=T) cat("Long-term effect of switch to SMD:",long.term,"% \n \n") outtable$OLS.cfe <- c(round(matrix(t(summary(OLS.cfe)$coefficients[c(1:8,nrow(summary(OLS.cfe)$coefficients)),c(1,2)])),3), # coefficients and SEs "",length(OLS.cfe$residuals), # number of observations round(summary(OLS.cfe)$r.squared,3),"", # R^2 of model round(long.term,1)) # long term effect of SMD ################################################################################################################# ################################################################################################################# ## ## Save table of results to file modeltable-dataset.csv ## #write.table(outtable,file=paste("modeltable-",whichdata,".csv",sep=""),sep=",",row.names=F) ################################################################################################################# ################################################################################################################# ################################################################################################################# ## Yearly real prices in SMD and PR democracies: 5-year moving average. par(mar=c(5,4.5,0.5,0.5),mfrow=c(1,1)) prices.pr <- matrix(0,ncol=1,nrow=29) prices.smd <- matrix(0,ncol=1,nrow=29) for (i in 1972:2000) { prices.pr[i-1971] <- mean(world$P[(world$year==i) & (world$SMD==0)],na.rm=T) prices.smd[i-1971] <- mean(world$P[(world$year==i) & (world$SMD==1)],na.rm=T) } ma.prices.pr <- matrix(0,ncol=1,nrow=25) ma.prices.smd <- matrix(0,ncol=1,nrow=25) for (i in 1:25) { ma.prices.pr[i] <- mean(prices.pr[i:(i+4)]) ma.prices.smd[i] <- mean(prices.smd[i:(i+4)]) } yvars <- cbind(ma.prices.pr,ma.prices.smd) matplot(names(table(world$year))[5:29],yvars,type="l",xlim=c(1970,2000),ylim=c(0,120), lty=c(2,1),col=1,cex.axis=1.3,cex.lab=1.3,las=1,lwd=c(2,2),xlab="Year",ylab="Real Price Index") text(1987,80,"PR",cex=1.3) text(1987,50,"SMD",cex=1.3) #savePlot(filename=paste(whichdata,"_ma_prices",sep=""),type="eps") ################################################################################################################# ## Yearly "baseline" prices in SMD and PR democracies, after controlling for GDP, changes ## in the exchange rate, population, and import levels: 5-year moving average. Vertical axis is the sum ## of GEE model intercept plus coefficient on electoral system-year dummy. cat("\n *************************************************** ") cat("\n *************** GEE Model: Year FEs *************** \n \n") year.smd <- matrix(0,nrow=nrow(world),ncol=(1+max(world$year)-min(world$year))) year.pr <- matrix(0,nrow=nrow(world),ncol=(1+max(world$year)-min(world$year))) for (i in (min(world$year)):(max(world$year))) { year.smd[((world$year==i) & (world$SMD==1)),(i-min(world$year)+1)] <- 1 year.pr[((world$year==i) & (world$SMD==0)),(i-min(world$year)+1)] <- 1 } year.smd <- data.frame(year.smd) year.pr <- data.frame(year.pr) names(year.smd) <- c("smd1972","smd1973","smd1974","smd1975","smd1976","smd1977","smd1978","smd1979","smd1980","smd1981","smd1982","smd1983","smd1984","smd1985","smd1986","smd1987","smd1988","smd1989","smd1990","smd1991","smd1992","smd1993","smd1994","smd1995","smd1996","smd1997","smd1998","smd1999","smd2000") names(year.pr) <- c("pr1972","pr1973","pr1974","pr1975","pr1976","pr1977","pr1978","pr1979","pr1980","pr1981","pr1982","pr1983","pr1984","pr1985","pr1986","pr1987","pr1988","pr1989","pr1990","pr1991","pr1992","pr1993","pr1994","pr1995","pr1996","pr1997","pr1998","pr1999","pr2000") world.yfe <- cbind(world,year.smd,year.pr) # run GEE exchangeable model with year fixed effects (YFE) GEE.yfe <- gee(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+SMD:ROOTGDPrc+ smd1974+smd1975+smd1976+smd1977+smd1978+smd1979+smd1980+smd1981+smd1982+smd1983+smd1984+smd1985+smd1986+smd1987+smd1988+smd1989+smd1990+smd1991+smd1992+smd1993+smd1994+smd1995+smd1996+smd1997+smd1998+smd1999+smd2000+ pr1974+pr1975+pr1976+pr1977+pr1978+pr1979+pr1980+pr1981+pr1982+pr1983+pr1984+pr1985+pr1986+pr1987+pr1988+pr1989+pr1990+pr1991+pr1992+pr1993+pr1994+pr1995+pr1996+pr1997+pr1998+pr1999+pr2000, id=cnum, data=world.yfe, corstr="exchangeable") cat("\n Estimates of over-time gap between real price levels in SMD/PR democracies: \n \n") baseline.smd <- coefficients(GEE.yfe)[2] + rbind(0,matrix(coefficients(GEE.yfe)[8:34])) baseline.pr <- rbind(0,matrix(coefficients(GEE.yfe)[35:61])) vcv <- GEE.yfe$robust.variance std.err <- sqrt(vcv[2,2]) # 1973 for (i in 8:34) { # 1974-2000 std.err <- c(std.err,sqrt(vcv[2,2]+vcv[i,i]+vcv[i+27,i+27]+2*(vcv[2,i]-vcv[2,i+27]-vcv[i,i+27]))) } baseline.diff <- data.frame(year=matrix(c(1973:2000)),SMD.minus.PR=baseline.smd-baseline.pr,std.err=std.err,z=(baseline.smd-baseline.pr)/std.err) baseline.diff$p <- round(2*(1-pnorm(abs(baseline.diff$z))),3) print(baseline.diff) cat("\n Average difference since 1973 of SMD-PR real price levels:", mean(baseline.diff[,2]), "\n \n") plot(baseline.diff$year,baseline.diff$SMD.minus.PR,ylim=c(-15,10),pch=19,cex=1.5,cex.axis=1.3,cex.lab=1.3,las=1,xlab="Year",ylab="Effect of SMD, 95% error bars") segments(baseline.diff$year,baseline.diff$SMD.minus.PR-1.96*baseline.diff$std.err,baseline.diff$year,baseline.diff$SMD.minus.PR+1.96*baseline.diff$std.err) segments(baseline.diff$year,baseline.diff$SMD.minus.PR-1.645*baseline.diff$std.err,baseline.diff$year,baseline.diff$SMD.minus.PR+1.645*baseline.diff$std.err,lwd=3) abline(0,0,lty=3) #savePlot(filename=paste(whichdata,"_baseline_prices_diff_rawSE",sep=""),type="eps") ma.baseline.diff <- matrix(NA,ncol=1,nrow=24) for (i in 1:24) { ma.baseline.diff[i] <- mean(baseline.diff[i:(i+4),2]) } plot(c(1977:2000),ma.baseline.diff,type="l",ylim=c(-8,3),xlim=c(1970,2000),lty=c(1,2),col=c(1,1),cex.axis=1.3,cex.lab=1.3,las=1,lwd=c(2,2),xlab="Year",ylab="Estimated SMD-PR Price Index Gap") points(c(1977:2000),ma.baseline.diff,pch=19,cex=1.5) abline(0,0,lty=3) #savePlot(filename=paste(whichdata,"_baseline_prices_diff_ma",sep=""),type="eps") cat("\n ************************************************************** \n \n ") ################################################################################################################# ## Do a Hausman specification test to check the assumption of the random effects model. cat("Hausman specification tests. \n \n") pdata.frame(world.nomiss,id="country",time="year",name="world.pd") pmods <- plm(P~SMD+LAGP+ROOTGDPrc+LNDXR+LNPOP+IMPORTS+USAINF+SMD:ROOTGDPrc,data=world.pd) print(summary(pmods)) pmods.nopop <- plm(P~SMD+LAGP+ROOTGDPrc+LNDXR+IMPORTS+USAINF+SMD:ROOTGDPrc,data=world.pd) print(summary(pmods.nopop)) cat("\n ************************************************************** \n \n ") ## ## List of democracy country-years and systems in data set ## listcy <- function(world) { cat("Democracy years and electoral systems. \n \n") tab <- table(world$country,world$year) for (i in 1:nrow(tab)) { outtxt <- NULL foundst <- foundend <- F typenext <- typest <- NA demyrs <- 0 for (j in 1:ncol(tab)) { if (tab[i,j]==1) { demyrs <- demyrs + 1 idx <- ((matrix(world$country)==rownames(tab)[i]) & (world$year==colnames(tab)[j])) } if (!foundst) { if (tab[i,j]==1) { styr <- colnames(tab)[j] foundst <- T typest <- ifelse(world$MIX[idx]==1,ifelse(world$SMD[idx]==1,"SMD (Mixed)","PR (Mixed)"),ifelse(world$SMD[idx]==1,"SMD","PR")) if (j!=ncol(tab)) { idxnext <- ((matrix(world$country)==rownames(tab)[i]) & (world$year==colnames(tab)[j+1])) typenext <- ifelse(world$MIX[idxnext]==1,ifelse(world$SMD[idxnext]==1,"SMD (Mixed)","PR (Mixed)"),ifelse(world$SMD[idxnext]==1,"SMD","PR")) if (typest!=typenext) { endyr <- colnames(tab)[j] foundend <- T } } } } else { if (j==ncol(tab)) { endyr <- "2000" foundend <- T } else { if (tab[i,j+1]==1) { idxnext <- ((matrix(world$country)==rownames(tab)[i]) & (world$year==colnames(tab)[j+1])) typenext <- ifelse(world$MIX[idxnext]==1,ifelse(world$SMD[idxnext]==1,"SMD (Mixed)","PR (Mixed)"),ifelse(world$SMD[idxnext]==1,"SMD","PR")) if (typest!=typenext) { endyr <- colnames(tab)[j] foundend <- T } } else { endyr <- colnames(tab)[j] foundend <- T } } } if (foundend) { if (styr==endyr) { outtxt <- paste(outtxt,typest," ",styr,", ",sep="") } else { outtxt <- paste(outtxt,typest," ",styr,"-",endyr,", ",sep="") } foundst <- foundend <- F } } cat(rownames(tab)[i]," (",demyrs,"): ",outtxt,"\n",sep="") } } #listcy(world) if (F) { par(mar=c(4,4,2,1),ask=T) changers <- c("France","Italy","Japan","New Zealand","Latvia","Mongolia","Poland","Sri Lanka","Thailand","Ukraine") for (i in 1:length(changers)) { cur <- world[world$country==changers[i],] plot(cur$year,cur$P,type="l",main=changers[i],xlab="Year",ylab="Real price level") for (j in 2:nrow(cur)) { if (cur$SMD[j] != cur$SMD[j-1]) { if (cur$SMD[j]==1) { abline(v=cur$year[j]) # regular = switching to SMD } else { abline(v=cur$year[j],lwd=3) # bold = switching to PR } } } } par(ask=F) } ## ## End of program ##