#graphs_compare.R
#
#This file produces the figures comparing the baseline and "no China" transition paths in
# the paper "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)
#note that the "_total" files are produced from the fm_transition and tf_transition programs,
#and must have this suffix manually added
load(file="graphout_fm_total.rda")
graphout.fm.total=graphout.fm
load(file="graphout_tf_total.rda")
graphout.tf.total=graphout.tf
load(file="graphout_tf.rda")
load(file="graphout_fm.rda")

T = dim(graphout.tf)[1]-1
assign('numyrs',10)

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

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

parameters.tot.imports=read.csv('outdata_total_imports.csv')
parameter.tot.imports=parameters.tot.imports[2:dim(parameters.tot.imports)[2]]

#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]))

x = 100*(exp( (1-beta) *(parameter.tot.imports$Wfm-Wfm)  )-1)
xstar = 100*(exp( (1-beta) *(parameter.tot.imports$Wstarfm-Wstarfm)  )-1)

xtf = 100*(exp( (1-beta) *(parameter.tot.imports$Wtf-Wtf)  )-1)
xstartf = 100*(exp( (1-beta) *(parameter.tot.imports$Wstartf-Wstartf)  )-1)


sink(file="welfare_comparative.txt",append=FALSE)
cat('  ','\n')
cat('Northern Welfare Gains from China, FM = ',x,'\n')
cat('Southern Welfare Gains from China, FM = ',xstar,'\n')

cat('Northern Welfare Gains from China, TF = ',xtf,'\n')
cat('Southern Welfare Gains from China, TF = ',xstartf,'\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)

china.col=3
china.sym=2
growth.lim = c(1.9,2.9)
xrate.lim= c(0.4,0.75)

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

plot(1:(T+1),100*((1+graphout.tf.total$g)^(1/numyrs)-1),type="l",main="Variety Growth",
     ylim=growth.lim,xlab="Period",ylab="Annual %",col=1,lwd=lwidnum)
points(1:(T+1),100*((1+graphout.tf.total$g)^(1/numyrs)-1),pch=1,col=1)
lines(1:(T+1),100*((1+graphout.tf$g)^(1/numyrs)-1),col=china.col,lwd=lwidnum)
points(1:(T+1),100*((1+graphout.tf$g)^(1/numyrs)-1),pch=china.sym,col=china.col)
lines(1:(T+1),100*((1+g.phi.prime.total)^(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)
legend("topright",c("Observed","No China Growth"),col=c(1,china.col),lty=1,lwd=lwidnum,bty="n",pch=c(1,china.sym))

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


dev.off()

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

plot(1:(T+1),100*((1+graphout.fm.total$g)^(1/numyrs)-1),type="l",main="Variety Growth",xlab="Period",ylab="Annual %",col=1,lwd=lwidnum)
points(1:(T+1),100*((1+graphout.fm.total$g)^(1/numyrs)-1),pch=1,col=1)
lines(1:(T+1),100*((1+graphout.fm$g)^(1/numyrs)-1),col=china.col,lwd=lwidnum)
points(1:(T+1),100*((1+graphout.fm$g)^(1/numyrs)-1),pch=china.sym,col=china.col)
lines(1:(T+1),100*((1+g.phi.prime.total)^(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)
legend("bottomright",c("Observed","No China Growth"),col=c(1,china.col),lty=1,lwd=lwidnum,bty="n",pch=c(1,china.sym))

plot(1:(T+1),graphout.fm.total$q,type="l",main="Exchange Rate",xlab="Period",ylab="North/South",col=1,lwd=lwidnum)
points(1:(T+1),graphout.fm.total$q,pch=1,col=1)
lines(1:(T+1),graphout.fm$q,col=china.col,lwd=lwidnum)
points(1:(T+1),graphout.fm$q,pch=china.sym,col=china.col)
lines(1:(T+1),q.phi.prime.total*ones((T+1),1),col=4,lwd=lwidnum)
lines(1:(T+1),q.phi*ones((T+1),1),col=2,lty=2,lwd=lwidnum)

dev.off()
