set scheme s1mono
insheet using aggregate.csv, clear
gen gdpGrow=log(gdp[_n+1]/gdp)
gen mfpGrow=log(mfp[_n+1]/mfp)
gen gdpPerCapGrow=log(gdp_real[_n+1]/gdp_real*population/population[_n+1])
save aggregate, replace

use gra, replace
drop if sic=="2911" | sic=="5172" | sic=="1311" |  sic=="1389" | substr(sic,1,3)=="490" | substr(sic,1,3)=="491" | substr(sic,1,3)=="492" | substr(sic,1,3)=="493" | sic=="4940" | substr(sic,1,1)=="6"
drop if loc~="USA"
drop if sale<0 | sale==. | emp<0 | emp==.

ren fyear year
gen sic3=substr(sic,1,2)  /* Get the three digit industry code */
destring sic3, replace force
merge n:1 year using aggregate
sort gvkey year
replace sale=sale/gdp*gdp_real  /* convert nominal sales to real sales */
gen lp=log(sale/emp)  /* labor productivity*/
gen lpGrow=lp[_n+1]-lp if year[_n+1]-year==1 & gvkey[_n+1]==gvkey 
gen salesGrow=log(sale[_n+1]/sale) if year[_n+1]-year==1 & gvkey[_n+1]==gvkey
gsort year -sale
by year: gen top=_n
replace top=(top<=100)   /* indicator, equals 1 if it is one of the hundred largest firms of that year  */
gen lpGrowTop=lpGrow if top==1
by year: egen mnProdGrow=mean(lpGrowTop)  
bys year sic3: egen mnProdGrow3=mean(lpGrowTop)

gen lpGrowResid=lpGrowTop-mnProdGrow
gen lpGrowResid3=lpGrowTop-mnProdGrow3

replace lpGrowResid=.2 if lpGrowResid>.2 & lpGrowResid~=.  /* censor labor productivity growth rates at 20% */ 
replace lpGrowResid=-.2 if lpGrowResid<-.2 & lpGrowResid~=.
replace lpGrowResid3=.2 if lpGrowResid3>.2 & lpGrowResid3~=.
replace lpGrowResid3=-.2 if lpGrowResid3<-.2 & lpGrowResid3~=.

gen weight=sale/gdp_real/1000
drop if _merge==2
drop _merge
replace year=year+1
bys year: egen granular=sum(lpGrowResid*weight)
bys year: egen granularInd=sum(lpGrowResid3*weight)

by year: gen fi=1 if _n==1  /* Keep one observation per year */
keep if fi==1
tset year
/* main regressions */
reg gdpPerCapGrow granular l.granular if year>=1952 & year<=2007  
reg gdpPerCapGrow granularInd l.granularInd if year>=1952 & year<=2007
gen combo=(_b[granularInd]*granularInd+_b[l.granularInd]*l.granularInd)
gen year2=mod(year,100)  /* get the last two digits of the year, so we can put it in the plot */
tostring year2, replace force
replace year2="0"+year2 if year>=2000 & year<=2009
graph twoway (scatter gdpPerCapGrow combo if year>=1952 & year<2008, mlabel(year2) mlabpos(0) mlabcolor(black) msymbol(i))  (scatter gdpPerCapGrow combo if year>=2008, mlabel(year2) mlabcolor(red) mlabpos(0) msymbol(i))  ,   xtitle("Granular Residual") ytitle("GDP Per Capita Growth") legend(off)
graph export "granular.eps", replace