#read.R
#
# This file performs the calculations of the patent ratios for non-OECD countries
# and China reported in 
# the paper "Trapped Factors and China's Impact on Growth,"
# Bloom, Nicholas, Romer, Paul, Terry, Stephen J., and Van Reenen, John (2013).
#
# The patent data was downloaded from the "NBER Patent Data Project,"
# at https://sites.google.com/site/patentdataproject/, covering 1976-2006
#
#If you have questions about the code, please contact Stephen Terry at sterry@stanford.edu.
#
#September 27, 2013


#read in NBER patent data, to examine OECD vs China patenting in the US
rm(list=ls())
library(matlab)
library(stats)
library(foreign)
options(digits=3)

#read in dta file with patent data from Brownwyn Hall's website
patents=read.dta(file="patsic06_mar09_ipc.dta")

#load patent assignee data
assign = read.dta("patassg.dta")
assign$patent = assign$patnum
assign$patnum=NULL

#merge the patent and patent assignee data
patents = merge(x=patents,y=assign,by=c("patent","pdpass"),all.x=TRUE,all.y=FALSE)
patents$cnt[patents$cnt==""]="US"
three.digits = nchar(patents$cnt)==3
three.digits[is.na(three.digits)]=FALSE
patents$cnt[three.digits]=substr(patents$cnt[three.digits],1,2)

#A PIECE OF CODE HERE THAT DOES THE SUBSAMPLING AND ONLY COUNTS A SINGLE PATENT ONCE
first.assignee = patents$assgnum==1
patents = patents[first.assignee,]

#now, read in formatted country codes and OECD membership data
country = read.csv(file="country_codes.csv")
country$cnt=country$code
country$code=NULL
country$country.name=country$name
country$name=NULL
country$oecd.country=country$oecd
country$oecd=NULL

#merge the country codes & OECD accession data with the patent info
patents=merge(x=patents,y=country[,c("country.name","cnt","oecd.country","oecd.year")],by="cnt",all.x=TRUE,all.y=FALSE)
sample=!is.na(patents$country.name)
patents=patents[sample,]

#now, for each patent determine if it was applied to from an OECD country (i.e. OECD identification comes from application date)
oecd.sample = patents$oecd.country==1
oecd.sample[is.na(oecd.sample)]=FALSE
patents$oecd.patent=NA
patents$oecd.patent[oecd.sample] = as.numeric(patents$appyear[oecd.sample]>=patents$oecd.year[oecd.sample])
patents$oecd.patent[is.na(patents$oecd.patent)]=0

#now, for each grant year, compute the fraction of patents from non-OECD countries, and from China
styr=1977
endyr=2006
numyrs=endyr-styr+1

patent.fractions = array(0,dim=c(numyrs,4))
patent.fractions = data.frame(patent.fractions)
colnames(patent.fractions)=c("appyear","NONOECD.FRACTION","CHINA.FRACTION","NON.US.FRACTION")

count=0
for (yrct in styr:endyr){
  count=count+1
  patent.fractions$appyear[count]=yrct
  yr.sample = patents$appyear==yrct
  numpatentsyr = sum(yr.sample,na.rm=TRUE)
  
  non.oecd.sample = yr.sample & (patents$oecd.patent==0)
  patent.fractions$NONOECD.FRACTION[count]=sum(non.oecd.sample,na.rm=TRUE)/numpatentsyr
  
  china.sample = yr.sample & (patents$cnt=="CN")
  patent.fractions$CHINA.FRACTION[count]=sum(china.sample,na.rm=TRUE)/numpatentsyr
  
  non.us.sample = yr.sample & (patents$cnt!="US")
  patent.fractions$NON.US.FRACTION[count]=sum(non.us.sample,na.rm=TRUE)/numpatentsyr
  
}

country.list=levels(factor(patents$cnt))
num.countries=length(country.list)

frac.patents.countries=array(0,dim=c(endyr-styr+1,num.countries))

pos=0
for (yrct in styr:endyr){
  pos = pos+1
  yr.sample = patents$appyear==yrct
  numpatentsyr = sum(yr.sample,na.rm=TRUE)
  cat(yrct,'\n')
  
  for (countryct in (1:num.countries)){
    countryval = country.list[countryct]  
    cat(countryval,' ')
    
    countryyrsamp = patents$cnt==countryval & yr.sample
    numcountryyr = sum(countryyrsamp,na.rm=TRUE)
    
    frac.patents.countries[pos,countryct]=numcountryyr/numpatentsyr
    
  }
    
}

frac.patents.countries=data.frame(frac.patents.countries)
colnames(frac.patents.countries)=country.list
frac.patents.countries=cbind((styr:endyr),frac.patents.countries)
colnames(frac.patents.countries)[1]="appyear"

write.csv(frac.patents.countries,file="frac_patents_countries.csv")

#compute 90's means (actually, the means from styr to endyr)
ninetiessmeans=c(0.0,0.0,0.0)
ninetiessample=patent.fractions$appyear>=styr & patent.fractions$appyear<=endyr
ninetiessmeans[1]=mean(100*patent.fractions$NON.US.FRACTION[ninetiessample],na.rm=TRUE)
ninetiessmeans[2]=mean(100*patent.fractions$NONOECD.FRACTION[ninetiessample],na.rm=TRUE)
ninetiessmeans[3]=mean(100*patent.fractions$CHINA.FRACTION[ninetiessample],na.rm=TRUE)

legendlist = list();
legendlist[[1]]=paste(paste("Total Foreign (mean = ",substr(as.character(ninetiessmeans[1]),1,4),sep=""),"%)",sep="")
legendlist[[2]]=paste(paste("Non-OECD (mean = ",substr(as.character(ninetiessmeans[2]),1,4),sep=""),"%)",sep="")
legendlist[[3]]=paste(paste("China (mean = ",substr(as.character(ninetiessmeans[3]),1,4),sep=""),"%)",sep="")

#plot the three different patent fractions
lwidnum=2.5
fsize=10
pdf("patent_fractions.pdf")
par(ps = fsize, cex = 1, cex.main = 1.5,cex.lab=1.5,cex.axis=1.25)
plot(patent.fractions$appyear,100*patent.fractions$NON.US.FRACTION,type="l",col=1,lty=1,lwd=lwidnum,main="US Patents from Foreign Countries, 1977-2006",
     ylab="Percent of all US Patents",ylim=c(0,100),xlab="Application Year")
points(patent.fractions$appyear,100*patent.fractions$NON.US.FRACTION,pch=1,col=1,)
lines(patent.fractions$appyear,100*patent.fractions$NONOECD.FRACTION,col=2,lty=1,lwd=lwidnum)
points(patent.fractions$appyear,100*patent.fractions$NONOECD.FRACTION,col=2,pch=2)
lines(patent.fractions$appyear,100*patent.fractions$CHINA.FRACTION,col=4,lty=1,lwd=lwidnum)
points(patent.fractions$appyear,100*patent.fractions$CHINA.FRACTION,col=4,pch=3)
legend("topright",cex=1.5,c(legendlist[[1]],legendlist[[2]],legendlist[[3]]),col=c(1,2,4),lty=c(1,1,1),pch=c(1,2,3),lwd=lwidnum,bty="n")
dev.off()

pdf("patent_fractions_China.pdf")
plot(patent.fractions$appyear,100*patent.fractions$CHINA.FRACTION,type="l",col=2,lwd=lwidnum,main="US Patents from China, 1977-2006",
     ylab="Percent of all US Patents",xlab="Application Year",ylim=c(0,1))
legend("topleft",c(paste(paste("Mean = ",substr(as.character(ninetiessmeans[3]),1,4),sep=""),"%",sep="")),col=2,lwd=lwidnum,bty="n")
dev.off()


save(file="patents.rda",list=ls())