#graphs.R
#
#This file produces the figures and performs welfare calculations in
#"Trapped Factors and China's Impact on Growth,"
#Bloom, Nicholas, Romer, Paul, Terry, Stephen J., and Van Reenen, John (2013).
#
#If you have questions about the code, please contact Stephen Terry at sterry@stanford.edu.
#
#September 27, 2013

rm(list=ls())
library(matlab)
options(digits=15)
load(file="graphout_tf.rda")
load(file="graphout_fm.rda")
T = dim(graphout.tf)[1]-1
assign('numyrs',10)


g.phi = graphout.fm$g[1]
g.phi.prime = graphout.fm$g[T+1]
q.phi=graphout.fm$q[1]
q.phi.prime=graphout.fm$q[T+1]
r.phi=graphout.fm$r[1]
r.phi.prime=graphout.fm$r[T+1]


scale = 1000/graphout.fm$M1[1]
numperbefore = 1
numperafter = 3

lwidnum=2

parameters = read.csv('parameters.csv')
assign('A1',graphout.fm$A[1])
assign('phi',parameters$phi)
assign('phi.prime',parameters$phi.prime)
assign('alpha',parameters$alpha)
assign('H',parameters$H)
assign('Hstar',parameters$Hstar)
assign('gamma',1/parameters$rho)
assign('eta',parameters$eta)
assign('IoverY.open.target',parameters$IoverY.open.target)
assign('IoverY.fully.mobile.target',parameters$IoverY.fully.mobile.target)
assign('beta',parameters$beta);  #discount rate


#get axis right
lengthgraphA = numperbefore + numperafter+1
graphA = zeros(lengthgraphA,5)
graphA[(numperbefore+1):lengthgraphA,1]=(1:(numperafter+1))
graphA[1:numperbefore,1] = linspace(-numperbefore+1,0,numperbefore)
graphA = graphA +1 

#FM patents
graphA[(numperbefore+1):lengthgraphA,2] = scale*graphout.fm$M1[1:(numperafter+1)]
for (t in (1:numperbefore)){
  graphA[numperbefore+1-t,2]=graphA[numperbefore+1-t+1,2]/(1+g.phi)
}

#industry 1 patents
graphA[(numperbefore+1):lengthgraphA,3] = scale*graphout.tf$M1[1:(numperafter+1)]
for (t in (1:numperbefore)){
  graphA[numperbefore+1-t,3]=graphA[numperbefore+1-t+1,3]/(1+g.phi)
}

#industry 2 patents
graphA[(numperbefore+1):lengthgraphA,4] = scale*graphout.tf$M2[1:(numperafter+1)]
for (t in (1:numperbefore)){
  graphA[numperbefore+1-t,4]=graphA[numperbefore+1-t+1,4]/(1+g.phi)
}



bar.data = c(scale*graphout.tf$M2[1],scale*graphout.tf$M1[2],0,scale*graphout.tf$M2[2])


pdf(file="micro_graph_text.pdf",pointsize=18)
barplot(bar.data[c(1,2,4)],border=NA,space = 0.5,names.arg=c("Pre-Shock","No Shock","Shocked"),
        ylab="New Patent Flows",xlab="Industry",col=c(1,4,4),ylim=c(800,2000),xpd=FALSE)
text(2.55,1600,"28.8%",font=2)
arrows(3.25,bar.data[2],3.25,bar.data[4],code=3,lwd=1.5,angle=90,length=.2)

dev.off()




#welfare calculations

#SET UP STORAGE
welfare.length = 600
u.mat=zeros(welfare.length,6)
betas=zeros(welfare.length,1)
u.mat[1:(T+1),1]=graphout.tf$u
u.mat[1:(T+1),2]=graphout.tf$ustar
u.mat[1:(T+1),3]=graphout.fm$u
u.mat[1:(T+1),4]=graphout.fm$ustar

#FIRST, COMPUTE THE WELFARE WITH NO SHOCK IN THE WEST AND THE EAST
Mprimeval = g.phi * A1; Mval = Mprimeval / (1+g.phi);
Rprimeval = (1-phi) * A1; Rval = Rprimeval / (1+g.phi);
Iprimeval = phi * A1; Ival = Iprimeval / (1+g.phi);

xMval = (1-alpha)^(2/alpha) * (1+r.phi)^(-1/alpha) * H
xMstarval = (1-alpha)^(2/alpha) * (1+r.phi)^(-1/alpha) * q.phi^(1/alpha) * Hstar

xRval = (1-alpha)^(1/alpha) * (1+r.phi)^(-1/alpha) * H
xRstarval = (1-alpha)^(1/alpha) * (1+r.phi)^(-1/alpha) * Hstar

xIval = (1-alpha)^(1/alpha) * (1+r.phi)^(-1/alpha) * q.phi^(-1/alpha) * H
xIstarval = (1-alpha)^(1/alpha) * (1+r.phi)^(-1/alpha) * Hstar

Zval = (g.phi^gamma)/(eta*(gamma-1)+1)

Yval = H^alpha
Yval = Yval * ( Mval * xMval^(1-alpha) + Rval * xRval^(1-alpha) + Ival * xIval^(1-alpha)  )

Cval = Yval - Mprimeval * (xMval + xMstarval) - Rprimeval * xRval - Zval

Ystarval = Hstar^alpha
Ystarval = Ystarval *  (Mval * xMstarval^(1-alpha) + Rval * xRstarval^(1-alpha) + Ival * xIstarval^(1-alpha)    )

Cstarval = Ystarval - Rprimeval * xRstarval - Iprimeval * (xIval + xIstarval)
#THEN, DISCOUNT ALL FLOWS

betas[1]=beta
u.mat[1,5]=Cval
u.mat[1,6]=Cstarval
for (t in (2:welfare.length)){
  
  u.mat[t,5:6] = u.mat[t-1,5:6]*(1+g.phi)
}

for (t in (2:welfare.length)){
  betas[t]=beta*betas[t-1]
}

for (t in ((T+2):welfare.length)) {
  
  u.mat[t,1:4]=(1+g.phi.prime)*u.mat[t-1,1:4]
  
}

Wtf = sum(betas*log(u.mat[,1]))
Wstartf = sum(betas*log(u.mat[,2]))
Wfm = sum(betas*log(u.mat[,3]))
Wstarfm = sum(betas*log(u.mat[,4]))
Wns = sum(betas*log(u.mat[,5]))
Wstarns = sum(betas*log(u.mat[,6]))

x = 100*(exp( (1-beta) *(Wtf-Wfm)  )-1)
xstar = 100*(exp( (1-beta) *(Wstartf-Wstarfm)  )-1)

xgains = 100*(exp( (1-beta) *(Wfm-Wns)  )-1)
xstargains=100*(exp( (1-beta) *(Wstarfm-Wstarns)  )-1)

xgainstf = 100*(exp( (1-beta) *(Wtf-Wns)  )-1)
xstargainstf=100*(exp( (1-beta) *(Wstartf-Wstarns)  )-1)

sink(file="welfare.txt",append=FALSE)
cat('  ','\n')
cat('Northern Cons Equivalent of trade in % = ',xgains,'\n')
cat('Southern Cons Equivalent of trade in % = ',xstargains,'\n')
cat('Northern Cons Equivalent of TF in % = ',x,'\n')
cat('Southern Cons Equivalent of TF in % = ',xstar,'\n')
cat('Northern Cons Equivalent of TF (total) in % = ',xgainstf,'\n')
cat('Southern Cons Equivalent of TF (total) in % = ',xstargainstf,'\n')
cat('  ','\n')

cat('Diff in Diff Patent Growth in % = ',100*log(graphout.tf$M2[2]/graphout.tf$M1[2]),'\n')
cat('  ','\n')

cat('Number of Years',numyrs,'\n')
cat('Pre-Shock phi = ',100*graphout.fm$phi[1],'\n')
cat('Post-Shock phi = ',100*graphout.fm$phi[2],'\n')
cat('Pre-Shock I/Yg = ',100*IoverY.open.target,'\n')
cat('Post-Shock I/Yg = ',100*IoverY.fully.mobile.target,'\n')
cat('Pre-Shock g (growth rate) = ',100*((1+g.phi)^(1/numyrs)-1),'\n')
cat('Post-Shock g (growth rate)= ',100*((1+g.phi.prime)^(1/numyrs)-1),'\n')
cat('LR Difference',100*((1+g.phi.prime)^(1/numyrs)-1)-100*((1+g.phi)^(1/numyrs)-1),'\n')
cat('Period 2 growth rate (fully mobile)',100*((1+graphout.fm$g[2])^(1/numyrs)-1),'\n')
cat('Period 2 growth rate (trapped factors)',100*((1+graphout.tf$g[2])^(1/numyrs)-1),'\n')
cat('SR Difference',100*((1+graphout.tf$g[2])^(1/numyrs)-1)-100*((1+graphout.fm$g[2])^(1/numyrs)-1),'\n')
cat('Pre-Shock q (x-rate) = ',q.phi,'\n')
cat('Post-Shock q (x-rate) = ',q.phi.prime,'\n')
cat('Pre-Shock r (interest rate)= ',100*((1+r.phi)^(1/numyrs)-1),'\n')
cat('Post-Shock r (interest rate) = ',100*((1+r.phi.prime)^(1/numyrs)-1),'\n')


sink(file=NULL)





#create the counterfactual decompositions for output in the shock period

xfun = function(p,alpha,H){
  obj = (1-alpha)^(1/alpha) * (p^(-1/alpha)) * H
}

#FM

pM2.fm = (1+graphout.fm$r[2])/(1-alpha)
pR2.fm = (1+graphout.fm$r[2])
pI2.fm = graphout.fm$q[2] * (1+graphout.fm$rstar[2])

xM2.fm = xfun(pM2.fm,alpha,H)
xR2.fm = xfun(pR2.fm,alpha,H)
xI2.fm = xfun(pI2.fm,alpha,H)

R2.fm = (1-graphout.fm$phi[2]) * graphout.fm$A[1]
I2.fm = graphout.fm$phi[2] * graphout.fm$A[1]
M2.fm = graphout.fm$M[2]

Y2.fm = M2.fm * (xM2.fm^(1-alpha)) + R2.fm * (xR2.fm^(1-alpha)) + I2.fm  * (xI2.fm^(1-alpha))
Y2.fm = (H^alpha) * Y2.fm


#no price effect counterfactual
R2noprice.fm = (1-graphout.fm$phi[1]) * graphout.fm$A[1]
I2noprice.fm = graphout.fm$phi[1] * graphout.fm$A[1]
M2noprice.fm = M2.fm

Y2noprice.fm = M2noprice.fm * (xM2.fm^(1-alpha)) + R2noprice.fm * (xR2.fm^(1-alpha)) + I2noprice.fm  * (xI2.fm^(1-alpha))
Y2noprice.fm = (H^alpha) * Y2noprice.fm

#no variety effect counterfactual
R2novariety.fm = R2noprice.fm
I2novariety.fm = I2noprice.fm
M2novariety.fm = graphout.fm$A[1] * graphout.fm$g[1]

Y2novariety.fm = M2novariety.fm * (xM2.fm^(1-alpha)) + R2novariety.fm * (xR2.fm^(1-alpha)) + I2novariety.fm  * (xI2.fm^(1-alpha))
Y2novariety.fm = (H^alpha) * Y2novariety.fm


PG.fm = 100*((Y2.fm/Y2noprice.fm)-1)
TG.fm = 100*((Y2.fm/graphout.fm$Y[1])-1)

cat('FM Case Counterfactuals','\n')
countervec.fm = c(Y2novariety.fm,Y2noprice.fm,Y2.fm)/Y2.fm
cat(countervec.fm,'\n')
cat('% due to R-->I price',100*(countervec.fm[3]-countervec.fm[2])/(countervec.fm[3]-countervec.fm[1]),'\n')
cat(' ','\n')

#TF

pM12.tf = tf.allocations$mu1[2] * (1+graphout.tf$r[2])/(1-alpha)
pM22.tf = tf.allocations$mu2[2] * (1+graphout.tf$r[2])/(1-alpha)
pR2.tf = (1+graphout.tf$r[2]); #note that mu1 and mu2 > 1/3 so monopoly price is bound by fast copiers
pI2.tf = graphout.tf$q[2] * (1+graphout.tf$rstar[2])

xM12.tf = xfun(pM12.tf,alpha,H)
xM22.tf = xfun(pM22.tf,alpha,H)
xR2.tf = xfun(pR2.tf,alpha,H)
xI2.tf = xfun(pI2.tf,alpha,H)

M12.tf = graphout.tf$M1[2]
M22.tf = graphout.tf$M2[2]
R2.tf = (1 - graphout.tf$phi[2] )* graphout.tf$A[1]
I2.tf = graphout.tf$phi[2] * graphout.tf$A[1]

Y2.tf = M12.tf * (xM12.tf^(1-alpha)) + M22.tf * (xM22.tf^(1-alpha)) + R2.tf * (xR2.tf^(1-alpha)) + I2.tf * (xI2.tf^(1-alpha))
Y2.tf = Y2.tf * (H^alpha)

#no price effect counterfactual

Mnoprice12.tf = M12.tf
Mnoprice22.tf = M22.tf
R2noprice.tf = (1-graphout.tf$phi[1]) * graphout.tf$A[1]
I2noprice.tf = graphout.tf$phi[1] * graphout.tf$A[1]



Y2noprice.tf = Mnoprice12.tf * (xM12.tf^(1-alpha)) + Mnoprice22.tf * (xM22.tf^(1-alpha)) +
  R2noprice.tf * (xR2.tf^(1-alpha)) + I2noprice.tf * (xI2.tf^(1-alpha))
Y2noprice.tf = Y2noprice.tf * (H^alpha)

#no variety effect counterfactual
R2novariety.tf = R2noprice.tf
I2novariety.tf = I2noprice.tf
Mnovariety12.tf = M12.tf * graphout.tf$g[1]/graphout.tf$g[2]
Mnovariety22.tf = M22.tf * graphout.tf$g[1]/graphout.tf$g[2]

Y2novariety.tf = Mnovariety12.tf * (xM12.tf^(1-alpha)) + Mnovariety22.tf * (xM22.tf^(1-alpha)) +
  R2novariety.tf * (xR2.tf^(1-alpha)) + I2novariety.tf * (xI2.tf^(1-alpha))
Y2novariety.tf = Y2novariety.tf * (H^alpha)

cat('TF Case Counterfactuals','\n')
countervec.tf= c(Y2novariety.tf,Y2noprice.tf,Y2.tf)/Y2.tf
cat(countervec.tf,'\n')
cat('% due to R-->I price',100*(countervec.tf[3]-countervec.tf[2])/(countervec.tf[3]-countervec.tf[1]),'\n')
cat(' ','\n')

#now, do the SS counterfactuals
rss = graphout.fm$r[1]
qss = graphout.fm$q[1]
gss = graphout.fm$g[1]
pMss = (1+rss)/(1-alpha)
pRss =  (1+rss)
pIss = qss * (1+rss)

xMss = xfun(pMss,alpha,H)
xRss = xfun(pRss,alpha,H)
xIss = xfun(pIss,alpha,H)

Ass = 1
Assmin1 = Ass/(1+gss)
Mss = Assmin1 * gss
Rss = Assmin1 * (1-phi)
Iss = Assmin1 * phi
Issmin1 = Iss/(1+gss)

#construct SS output
Yss = Mss * xMss^(1-alpha)  + Rss * xRss^(1-alpha) + Iss * xIss^(1-alpha)
Yss = H^alpha * Yss

#construct the SS output in the no R->I conversion case, holding total number of varieties constant
Mssnoprice=Mss
Rssnoprice=Assmin1 - Issmin1
Issnoprice=Issmin1

Yssnoprice = Mssnoprice * xMss^(1-alpha)  + Rssnoprice * xRss^(1-alpha) + Issnoprice * xIss^(1-alpha)
Yssnoprice = Yssnoprice * H^alpha

#construct no new varieties output
Yssnovariety = Yssnoprice - (H^alpha) * xMss^(1-alpha) * Mss

cat('SS Case Counterfactual','\n')
countervec.ss = c(Yssnovariety,Yssnoprice,Yss)/Yss
cat(countervec.ss,'\n')
cat('% due to R-->I price',100*(countervec.ss[3]-countervec.ss[2])/(countervec.ss[3]-countervec.ss[1]),'\n')
cat(' ','\n')






growth.lim = c(1.9,2.9)
xrate.lim= c(0.4,0.75)

pdf(file="tf_trans4.pdf")
par(mfrow=c(2,2))


plot(1:(T+1),100*((1+graphout.tf$g)^(1/numyrs)-1),type="l",main="Variety Growth",xlab="Period",ylab="Annual %",col=1,lwd=lwidnum,ylim=growth.lim)
points(1:(T+1),100*((1+graphout.tf$g)^(1/numyrs)-1),pch=1)
lines(1:(T+1),100*((1+g.phi.prime)^(1/numyrs)-1)*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),100*((1+g.phi)^(1/numyrs)-1)*ones((T+1),1),col=2,lty=2,lwd=lwidnum)

plot(1:(T+1),graphout.tf$q,type="l",main="Southern Terms of Trade",xlab="Period",ylab="North/South",col=1,lwd=lwidnum,ylim=xrate.lim)
points(1:(T+1),graphout.tf$q,pch=1)
lines(1:(T+1),q.phi.prime*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),q.phi*ones((T+1),1),col=2,lty=2,lwd=lwidnum)

plot(1:(T+1),c(100*((1+g.phi)^(1/numyrs)-1),100*((graphout.tf$Y[2:(T+1)]/graphout.tf$Y[1:T])^(1/numyrs)-1)),type="l",main="Northern Output Growth",xlab="Period",
     ylab="Annual %",col=1,lwd=lwidnum,ylim=growth.lim)
points(1:(T+1),c(100*((1+g.phi)^(1/numyrs)-1),100*((graphout.tf$Y[2:(T+1)]/graphout.tf$Y[1:T])^(1/numyrs)-1)),pch=1)
lines(1:(T+1),100*((1+g.phi.prime)^(1/numyrs)-1)*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),100*((1+g.phi)^(1/numyrs)-1)*ones((T+1),1),col=2,lty=2,lwd=lwidnum)

plot(1:(T+1),c(100*((1+g.phi)^(1/numyrs)-1),100*((graphout.tf$Ystar[2:(T+1)]/graphout.tf$Ystar[1:T])^(1/numyrs)-1)),type="l",main="Southern Output Growth",xlab="Period",ylab="Annual %",
     col=1,lwd=lwidnum,ylim=growth.lim)
points(1:(T+1),c(100*((1+g.phi)^(1/numyrs)-1),100*((graphout.tf$Ystar[2:(T+1)]/graphout.tf$Ystar[1:T])^(1/numyrs)-1)),pch=1)
lines(1:(T+1),100*((1+g.phi.prime)^(1/numyrs)-1)*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),100*((1+g.phi)^(1/numyrs)-1)*ones((T+1),1),col=2,lty=2,lwd=lwidnum)


dev.off()


pdf(file="fm_trans4.pdf")
par(mfrow=c(2,2))

plot(1:(T+1),100*((1+graphout.fm$g)^(1/numyrs)-1),type="l",main="Variety Growth",xlab="Period",ylab="Annual %",col=1,
     lwd=lwidnum,ylim=growth.lim)
points(1:(T+1),100*((1+graphout.fm$g)^(1/numyrs)-1),pch=1)
lines(1:(T+1),100*((1+g.phi.prime)^(1/numyrs)-1)*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),100*((1+g.phi)^(1/numyrs)-1)*ones((T+1),1),col=2,lty=2,lwd=lwidnum)

plot(1:(T+1),graphout.fm$q,type="l",main="Southern Terms of Trade",xlab="Period",ylab="North/South",col=1,lwd=lwidnum,
     ylim=xrate.lim)
points(1:(T+1),graphout.fm$q,pch=1)
lines(1:(T+1),q.phi.prime*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),q.phi*ones((T+1),1),col=2,lty=2,lwd=lwidnum)

plot(1:(T+1),c(100*((1+g.phi)^(1/numyrs)-1),100*((graphout.fm$Y[2:(T+1)]/graphout.fm$Y[1:T])^(1/numyrs)-1)),type="l",main="Northern Output Growth",xlab="Period",ylab="Annual %",col=1,
     lwd=lwidnum,ylim=growth.lim)
points(1:(T+1),c(100*((1+g.phi)^(1/numyrs)-1),100*((graphout.fm$Y[2:(T+1)]/graphout.fm$Y[1:T])^(1/numyrs)-1)),pch=1)
lines(1:(T+1),100*((1+g.phi.prime)^(1/numyrs)-1)*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),100*((1+g.phi)^(1/numyrs)-1)*ones((T+1),1),col=2,lty=2,lwd=lwidnum)

plot(1:(T+1),c(100*((1+g.phi)^(1/numyrs)-1),100*((graphout.fm$Ystar[2:(T+1)]/graphout.fm$Ystar[1:T])^(1/numyrs)-1)),type="l",main="Southern Output Growth",xlab="Period",ylab="Annual %",col=1,
     lwd=lwidnum,ylim=growth.lim)
points(1:(T+1),c(100*((1+g.phi)^(1/numyrs)-1),100*((graphout.fm$Ystar[2:(T+1)]/graphout.fm$Ystar[1:T])^(1/numyrs)-1)),pch=1)
lines(1:(T+1),100*((1+g.phi.prime)^(1/numyrs)-1)*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),100*((1+g.phi)^(1/numyrs)-1)*ones((T+1),1),col=2,lty=2,lwd=lwidnum)
dev.off()

