#PWT_PROCESS.R
#
# This file performs the calculation of the China --> OECD and non-OECD imports to GDP ratios
# used in calibration of the paper
# the paper "Trapped Factors and China's Impact on Growth,"
# Bloom, Nicholas, Romer, Paul, Terry, Stephen J., and Van Reenen, John (2013).
#
# The PWT data was downloaded on April 2, 2013
#
#If you have questions about the code, please contact Stephen Terry at sterry@stanford.edu.
#
#September 27, 2013

rm(list=ls())
library(foreign)
library(matlab)
library(stats)
options(digits=2)
#read in the PWT data
RAW.PWT=read.csv(file="pwt71_data.csv")

sample=!is.na(RAW.PWT$ppp)&!is.na(RAW.PWT$tcgdp)&!is.na(RAW.PWT$XRAT)
PWT=RAW.PWT[sample,]

#read in the ISO codes and OECD identifiers
ISO.CODES=read.csv(file="iso_codes.csv")
ISO.CODES$isocode=ISO.CODES$ISO3; ISO.CODES$ISO3=NULL; ISO.CODES$ISO2=NULL

#first, merge the ISO codes and the PWT data
PWT=merge(x=PWT,y=ISO.CODES,by="isocode",all.x=TRUE,all.y=FALSE)
country.year.order = order(PWT$isocode,PWT$year,decreasing=FALSE)
PWT=PWT[country.year.order,]
#then, provide an indicator for whether a particular data year is an OECD year or not
oecd.sample=PWT$OECD.COUNTRY==1
oecd.sample[is.na(oecd.sample)]=FALSE
PWT$OECD.OBS=0
PWT$OECD.OBS[oecd.sample]=PWT$year[oecd.sample]>=PWT$OECD.YEAR[oecd.sample]

#now, compute the ratio of ppp exchange rate to mkt exchange rate (if a dollar buys 3 RMB under PPP, but 8 at market, then 3/8 is the ratio)
PWT$ppp.ratio = PWT$ppp/PWT$XRAT
PWT$nominal.GDP = PWT$ppp.ratio*PWT$tcgdp

#now, for each year compute OECD nominal GDP, non-OECD nominal GDP, China nominal GDP, US nominal GDP
styr=1950
endyr=2010

numyrs=endyr-styr+1

SUMS.DATA=array(0,dim=c(numyrs,8))
SUMS.DATA=data.frame(SUMS.DATA)
colnames(SUMS.DATA)=c("year","OECD.NGDP","NON.OECD.NGDP","CHINA.NGDP","US.NGDP","OECD.PPP.RAT","NON.OECD.PPP.RAT","CHINA.PPP.RAT")

count=0
for (yrct in styr:endyr){
  count=count+1
  SUMS.DATA$year[count]=yrct
  
  yr.sample=PWT$year==yrct
  
  oecd.yr.sample=oecd.sample&yr.sample
  SUMS.DATA$OECD.NGDP[count]=sum(PWT$nominal.GDP[oecd.yr.sample],na.rm=TRUE)
  SUMS.DATA$OECD.PPP.RAT[count]=sum(PWT$nominal.GDP[oecd.yr.sample] * PWT$ppp.ratio[oecd.yr.sample]) / SUMS.DATA$OECD.NGDP[count]
  
  
  non.oecd.yr.sample = (!oecd.sample)&yr.sample
  SUMS.DATA$NON.OECD.NGDP[count]=sum(PWT$nominal.GDP[non.oecd.yr.sample],na.rm=TRUE)
  SUMS.DATA$NON.OECD.PPP.RAT[count]=sum(PWT$nominal.GDP[non.oecd.yr.sample] * PWT$ppp.ratio[non.oecd.yr.sample]) / SUMS.DATA$NON.OECD.NGDP[count]
  
  
  china.yr.sample = (PWT$isocode=="CHN")&yr.sample
  SUMS.DATA$CHINA.NGDP[count]=sum(PWT$nominal.GDP[china.yr.sample],na.rm=TRUE)
  SUMS.DATA$CHINA.PPP.RAT[count]=sum(PWT$nominal.GDP[china.yr.sample] * PWT$ppp.ratio[china.yr.sample]) / SUMS.DATA$CHINA.NGDP[count]
  
  us.yr.sample = (PWT$isocode=="USA")&yr.sample
  SUMS.DATA$US.NGDP[count]=sum(PWT$nominal.GDP[us.yr.sample],na.rm=TRUE)
  
}

SUMS.DATA$WORLD.NGDP=SUMS.DATA$OECD.NGDP+SUMS.DATA$NON.OECD.NGDP


#now, load in OECD data, with OECD gross output and non-OECD --> OECD imports
units.adjust=1000000
OECD.DATA=read.csv(file="non_oecd_imports.csv")
OECD.DATA$NON.OECD.IMPORTS=OECD.DATA$NON.OECD.IMPORTS/units.adjust
OECD.DATA$OECD.GROSS.OUTPUT=OECD.DATA$OECD.GROSS.OUTPUT/units.adjust

SUMS.DATA=merge(x=SUMS.DATA,y=OECD.DATA,by="year",all.x=TRUE,all.y=FALSE)

#now, compute non-OECD imports to OECD gross and value added output ratios
SUMS.DATA$OECD.IMPORTS.TO.GROSS=SUMS.DATA$NON.OECD.IMPORTS/SUMS.DATA$OECD.GROSS.OUTPUT
SUMS.DATA$OECD.IMPORTS.TO.GDP=SUMS.DATA$NON.OECD.IMPORTS/SUMS.DATA$OECD.NGDP
SUMS.DATA$post.shock = as.numeric(SUMS.DATA$year>=2002)

pre.shock.gross = mean(SUMS.DATA$OECD.IMPORTS.TO.GROSS[SUMS.DATA$post.shock==0],na.rm=TRUE)
post.shock.gross = mean(SUMS.DATA$OECD.IMPORTS.TO.GROSS[SUMS.DATA$post.shock==1],na.rm=TRUE)
pre.shock.gdp = mean(SUMS.DATA$OECD.IMPORTS.TO.GDP[SUMS.DATA$post.shock==0],na.rm=TRUE)
post.shock.gdp = mean(SUMS.DATA$OECD.IMPORTS.TO.GDP[SUMS.DATA$post.shock==1],na.rm=TRUE)


#post 1992 Chinese reforms average of PPP x-rate
post.reform=SUMS.DATA$year>=1992
post.reform.avg = mean(SUMS.DATA$CHINA.PPP.RAT[post.reform]/SUMS.DATA$OECD.PPP.RAT[post.reform],na.rm=TRUE)

#create graph of Chinese and non-OECD PPP x-rates to OECD values

SUMS.DATA$NON.OECD.TO.OECD.PPP.RAT.ADJUSTED = SUMS.DATA$NON.OECD.PPP.RAT/SUMS.DATA$OECD.PPP.RAT
SUMS.DATA$NON.OECD.TO.OECD.PPP.RAT.ADJUSTED[SUMS.DATA$year>=1990 & SUMS.DATA$year<=1994] = NA
non.oecd.post.reform.avg = mean(SUMS.DATA$NON.OECD.TO.OECD.PPP.RAT.ADJUSTED[post.reform],na.rm=TRUE)


lwidnum=2.5
pdf(file="CHINESE_MKT_TO_PPP.pdf")
plot(SUMS.DATA$year,SUMS.DATA$CHINA.PPP.RAT,ylim=c(0,1.0),type="l",col=4,lwd=lwidnum,xlim=c(1990,2010),
     main="Chinese to OECD Market to PPP Adjustment",ylab="Mkt/PPP x-rate, relative to OECD Mkt/PPP x-rate",xlab="Year")
text(2000,1.75,paste("Post-1992 Average = ",substr(as.character(post.reform.avg),1,4),sep=""))
dev.off()

legendlist=list()
legendlist[[1]]=paste("China, ",paste("Post-1992 Average = ",substr(as.character(post.reform.avg),1,4),sep=""),sep="")
legendlist[[2]]=paste("Non-OECD, ",paste("Post-1992 Average = ",substr(as.character(non.oecd.post.reform.avg),1,4),sep=""),sep="")
pdf(file="NON_OECD_MKT_TO_PPP.pdf")
plot(SUMS.DATA$year,SUMS.DATA$CHINA.PPP.RAT,ylim=c(0,1.5),type="l",col=4,lwd=lwidnum,xlim=c(1975,2010),
     main="Market to PPP Adjustment",ylab="Mkt/PPP x-rate, relative to OECD Mkt/PPP x-rate",xlab="Year")
lines(SUMS.DATA$year,SUMS.DATA$NON.OECD.TO.OECD.PPP.RAT.ADJUSTED,col=3,lwd=lwidnum)
legend("topright",c(legendlist[[1]],legendlist[[2]]),col=c(4,3),lwd=lwidnum,bty="n")
dev.off()

#now, load imports solely from China to the OECD
china.imports = read.csv(file="china_imports.csv")
china.imports=merge(x=china.imports,y=ISO.CODES,by="isocode",all.x=TRUE,all.y=FALSE)

imp.styr=min(china.imports$year)
imp.endyr=max(china.imports$year)

SUMS.DATA$china.imports=NA

for (yrct in (imp.styr:imp.endyr)){
  yr.samp = (china.imports$year==yrct)&(china.imports$OECD.YEAR<=china.imports$year)
  
  sum.dum = sum(china.imports$China.Imports[yr.samp],na.rm=TRUE)
  
  sums.data.obs=SUMS.DATA$year==yrct
    
  SUMS.DATA$china.imports[sums.data.obs]=sum.dum/1000; #need to convert OECD data ($K's USD) to PWT data ($mill USD)

}


pdf(file="IMPORT_RATIOS.pdf")
plot.data=SUMS.DATA[SUMS.DATA$year>=1997&SUMS.DATA$year<=2006,]
plot(plot.data$year,100*plot.data$NON.OECD.IMPORTS/plot.data$OECD.NGDP,col=1,type="l",
     ylim=c(0,7),lwd=lwidnum,main='Low-Cost Imports in the OECD',ylab="Percent of Total OECD GDP",xlab="Year")
lines(plot.data$year,100*plot.data$china.imports/plot.data$OECD.NGDP,col=4,lwd=lwidnum)
legend("topleft",c("Non-OECD","China"),col=c(1,4),lty=1,bty="n",lwd=lwidnum)
dev.off()

#output all of the average and ratios to a text file
sink(file="PWT.txt",append=FALSE)
cat('pre-shock imports to gross = ',pre.shock.gross,'\n')
cat('post-shock imports to gross = ',post.shock.gross,'\n')
cat('pre-shock imports to gdp = ',pre.shock.gdp,'\n')
cat('post-shock imports to gdp = ',post.shock.gdp,'\n')
cat('post-reform Chinese PPP-Adjustment Average = ',post.reform.avg,'\n')
cat('post-reform non-OECD PPP-Adjustment Average = ',non.oecd.post.reform.avg,'\n')
cat('','\n')
cat('pre-shock non-OECD imports to gdp, 1997 = ',SUMS.DATA$NON.OECD.IMPORTS[SUMS.DATA$year==1997]/SUMS.DATA$OECD.NGDP[SUMS.DATA$year==1997],'\n')
cat('pre-shock china imports to gdp, 1997 = ',SUMS.DATA$china.imports[SUMS.DATA$year==1997]/SUMS.DATA$OECD.NGDP[SUMS.DATA$year==1997],'\n')
cat('difference = ',SUMS.DATA$NON.OECD.IMPORTS[SUMS.DATA$year==1997]/SUMS.DATA$OECD.NGDP[SUMS.DATA$year==1997]-SUMS.DATA$china.imports[SUMS.DATA$year==1997]/SUMS.DATA$OECD.NGDP[SUMS.DATA$year==1997],'\n')
cat('post-shock non-OECD imports to gdp, 2006 = ',SUMS.DATA$NON.OECD.IMPORTS[SUMS.DATA$year==2006]/SUMS.DATA$OECD.NGDP[SUMS.DATA$year==2006],'\n')
cat('post-shock china imports to gdp, 2006 = ',SUMS.DATA$china.imports[SUMS.DATA$year==2006]/SUMS.DATA$OECD.NGDP[SUMS.DATA$year==2006],'\n')
cat('difference = ',SUMS.DATA$NON.OECD.IMPORTS[SUMS.DATA$year==2006]/SUMS.DATA$OECD.NGDP[SUMS.DATA$year==2006]-SUMS.DATA$china.imports[SUMS.DATA$year==2006]/SUMS.DATA$OECD.NGDP[SUMS.DATA$year==2006],'\n')
sink(file=NULL)
sink(file=NULL)
