#semiendog.R
#
#This file performs the calibration, computes the transition path, and 
#produces the figures for the semiendogenous growth robustness check version of the model 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(pso)
library(stats)
library(matlab)
library(compiler)
set.seed(1)
setCompilerOptions(optimize=3)

doplots=1

start.time=Sys.time()

#input parameters
numyrs=10
sigma=1
alpha=2/3
rho=0.5
eta=1
human.cap.ratio=2.96
beta=(1/1.02)^numyrs
N=2

#transition path parameters and normalization
T=100
A1=1
errortol=10^(-10)

#calibration targets and data
pop.growth = 0.0078
per.cap.growth = 0.0237
IoverY.preshock = 0.039
IoverY.postshock = 0.07

#determine the value of delta, as well as BGP growth rate
gH = (1+pop.growth)^numyrs-1
gA = (1+per.cap.growth)^numyrs-1
delta = 1 - rho*((log(1+gH))/(log(1+gA)))
rBGP = (1/beta)*(1+gH)*(1+gA)^sigma-1

#now, calibrate the values of phi and H1
import.targets = function(IN,paramvec,evaldum,flag.dum){
  
  if (flag.dum==1){ #this means that H1 and phi have already been calibrated, you're doing phi.prime
    
    phi = IN[1]
      
    sigma = paramvec[1]; alpha = paramvec[2];
    rho = paramvec[3]; eta = paramvec[4];
    human.cap.ratio = paramvec[5];  beta = paramvec[6]
    A1 = paramvec[7]; H1 = paramvec[8]
    Hstar1=paramvec[9]; gH = paramvec[10]
    gA = paramvec[11];  delta= paramvec[12]
    IoverY.target = paramvec[13]; 
    
    #constants
    PSI = (1-alpha)^((alpha-1)/(2-alpha))
    
    #BGP interest rates
    r = (1/beta)*(1+gH)*(1+gA)^sigma-1
    rstar = (1/beta)*(1+gH)*(1+gA)^sigma-1
    
    #BGP exchange rate
    q = ((phi/gA)*(1/human.cap.ratio))
    q = q^(alpha/(2-alpha)) * PSI
    
    #intensive margins
    xM = (1-alpha)^(2/alpha) * (1+r)^(-1/alpha)* H1
    xR = (1-alpha)^(1/alpha) * (1+r)^(-1/alpha)* H1
    xI = (1-alpha)^(1/alpha) * (1+rstar)^(-1/alpha)* q^(-1/alpha) * H1
    
    #extensive margins
    I = phi * A1/(1+gA)
    R = (1-phi)*A1/(1+gA)
    M = gA*A1/(1+gA)
    
    #create Northern output
    Y = H1^alpha
    Y = Y * ( M * xM^(1-alpha) + R * xR^(1-alpha) + I*xI^(1-alpha) )
    
    #create Northern imports
    pI = q * (1+rstar)
    Imports = I * pI * xI
    
    #create calibration moment
    IoverY=Imports/Y
    
    
    obj = (IoverY - IoverY.target)^2
    cat('error = ',obj,'\n')
    
    
    if (evaldum==0){ return(obj)}
    if (evaldum==1){ return(q)}
    
  }
  
  if (flag.dum==0){ #this means that phi and H1 have not yet been calibrated, you're doing both
    
    phi = IN[1]
    H1 = IN[2]

    sigma = paramvec[1]; alpha = paramvec[2];
    rho = paramvec[3]; eta = paramvec[4];
    human.cap.ratio = paramvec[5];  beta = paramvec[6]
    A1 = paramvec[7]; gH = paramvec[8]
    gA = paramvec[9];  delta= paramvec[10]
    IoverY.target = paramvec[11]; 
    
    
    #constants
    PSI = (1-alpha)^((alpha-1)/(2-alpha))
    OMEGA = alpha* ( 1 - alpha ) ^ ( ( 2 - alpha ) / alpha )
    gamma=1/rho
    
    #BGP interest rates
    r = (1/beta)*(1+gH)*(1+gA)^sigma-1
    rstar = (1/beta)*(1+gH)*(1+gA)^sigma-1
    
    #BGP exchange rate
    q = ((phi/gA)*(1/human.cap.ratio))
    q = q^(alpha/(2-alpha)) * PSI
    
    #intensive margins
    xM = (1-alpha)^(2/alpha) * (1+r)^(-1/alpha)* H1
    xR = (1-alpha)^(1/alpha) * (1+r)^(-1/alpha)* H1
    xI = (1-alpha)^(1/alpha) * (1+rstar)^(-1/alpha)* q^(-1/alpha) * H1
    
    #extensive margins
    I = phi * A1/(1+gA)
    R = (1-phi)*A1/(1+gA)
    M = gA*A1/(1+gA)
    
    #create Northern output
    Y = H1^alpha
    Y = Y * ( M * xM^(1-alpha) + R * xR^(1-alpha) + I*xI^(1-alpha) )
    
    #create Northern imports
    pI = q * (1+rstar)
    Imports = I * pI * xI
    
    #create calibration moment
    IoverY=Imports/Y
    
    #innov FOC
    LHSINNOV = gA^(gamma-1) * (1+gA)^((delta-1)/rho) * A1^((1-delta)/rho)
    RHSINNOV = OMEGA * (1+r)^(-1/alpha) * (H1 + q^(1/alpha)*human.cap.ratio*H1)
    
    
    obj = ((IoverY - IoverY.target)/IoverY)^2 + ((LHSINNOV-RHSINNOV)/RHSINNOV)^2
    cat('error = ',obj,'\n')
    
    
    if (evaldum==0){ return(obj)}
    if (evaldum==1){ return(q)}
    
  }
  
  
  
}


paramvec.preshock = c(sigma,alpha,rho,eta,
                      human.cap.ratio,beta,
                      A1,gH,gA,delta,
                      IoverY.preshock)

lb.preshock = c(10^-5,10^-5)
ub.preshock = c(0.3,5)

c.import.targets  = cmpfun(import.targets)

phi.list = psoptim(c(.075,.5),c.import.targets,lower=lb.preshock,upper=ub.preshock,paramvec=paramvec.preshock,evaldum=0,flag.dum=0)
phi = phi.list$par[1]
H1 = phi.list$par[2]
Hstar1=human.cap.ratio*H1
q.phi = import.targets(c(phi,H1),paramvec.preshock,evaldum=1,flag.dum=0)

paramvec.postshock = c(sigma,alpha,rho,eta,
                       human.cap.ratio,beta,
                       A1,H1,Hstar1,gH,gA,delta,
                       IoverY.postshock)

phi.prime.list = psoptim(0.075,c.import.targets,lower=10^-5,upper=0.3,paramvec=paramvec.postshock,evaldum=0,flag.dum=1)
phi.prime = phi.prime.list$par
q.phi.prime = c.import.targets(phi.prime,paramvec.postshock,evaldum=1,flag.dum=1)


fn.transition = function(pricevec,paramvec,evaldum){
  
  sigma = paramvec[1]; alpha = paramvec[2];
  rho = paramvec[3]; eta = paramvec[4];
  human.cap.ratio = paramvec[5];  beta = paramvec[6]
  A1 = paramvec[7]; H1 = paramvec[8]
  Hstar1=paramvec[9]; gH = paramvec[10]
  gABGP = paramvec[11];  delta= paramvec[12]
  T = paramvec[13]; rBGP = paramvec[14];
  q.phi = paramvec[15]; q.phi.prime=paramvec[16]
  phi = paramvec[17]; phi.prime = paramvec[18]
  
  #constants
  PSI = (1-alpha)^((alpha-1)/(2-alpha))
  OMEGA = alpha* ( 1 - alpha ) ^ ( ( 2 - alpha ) / alpha )
  gamma = 1/rho
  
  #create phi.matrix
  phivec = matrix(0,T+1,1)
  phivec[1]=phi
  phivec[2:(T+1)]=phi.prime
  
  #read in exchange rates
  q=matrix(0,T+1,1)
  q[1]=q.phi; q[T+1]=q.phi.prime
  q[2:T] = pricevec[1:(T-1)]
  
  #read in interest rates
  r = matrix(0,T+1,1)
  r[1]=rBGP;  r[T+1]=rBGP
  r[2:T]=pricevec[T:(2*(T-1))]
  
  rstar = matrix(0,T+1,1)
  rstar[1]=rBGP;  rstar[T+1]=rBGP
  rstar[2:T]=pricevec[(2*(T-1)+1):(3*(T-1))]
  
  #create the vector of human capital stocks
  H=matrix(0,T+1,1); H[1]=H1
  Hstar=matrix(0,T+1,1); Hstar[1]=Hstar1
  
  for (t in (2:(T+1))){
    H[t]=H[t-1]*(1+gH)
    Hstar[t]=Hstar[t-1]*(1+gH)
  }
  
  
  #now, simultaneously create the gA and A matrices
  gA = matrix(0,T+1,1); A=gA;
  gA[1]=gABGP; A[1]=A1;
  
  for (t in (2:(T))){
    
    gA[t] = OMEGA * (1+r[t])^(-1/alpha) * (H[t] + q[t]^(1/alpha) * Hstar[t]) * A[t-1]^((delta-1)/rho)
    gA[t] = gA[t]^(1/(gamma-1))
    
    A[t] = (1+gA[t])*A[t-1]
  }
  gA[T+1]=gABGP; A[T+1]=A[T]*(1+gA[T+1])
  
  #create the M, I, and R matrices
  M = matrix(0,T+1,1); I = M; R = M
  M[1] = gABGP * A1 /(1+gABGP)
  I[1] = phivec[1]*A1/(1+gABGP)
  R[1] = (1-phivec[1])*A1/(1+gABGP)
  
  for (t in (2:(T+1))){
    M[t] = gA[t]*A[t-1]
    I[t] = phivec[t]*A[t-1]
    R[t] = (1-phivec[t])*A[t-1]
  }
  
  #now, create intensive margins
  xM = (1-alpha)^(2/alpha) * (1+r)^(-1/alpha) * H
  xMstar = (1-alpha)^(2/alpha) * (1+r)^(-1/alpha) * q^(1/alpha) * Hstar
  
  xR = (1-alpha)^(1/alpha) * (1+r)^(-1/alpha) * H
  xRstar = (1-alpha)^(1/alpha) * (1+rstar)^(-1/alpha) * Hstar
  
  xI = (1-alpha)^(1/alpha) * (1+rstar)^(-1/alpha) * q^(-1/alpha) * H
  xIstar = (1-alpha)^(1/alpha) * (1+rstar)^(-1/alpha) * Hstar
  
  #now create output vectors for the North and the South
  Y = matrix(0,T+1,1); Ystar=Y;
  
  Y = H^alpha
  Y = Y * (M*xM^(1-alpha) + R*xR^(1-alpha) + I*xI^(1-alpha))
  
  Ystar = Hstar^alpha
  Ystar = Ystar *(M*xMstar^(1-alpha) + R*xRstar^(1-alpha) + I*xIstar^(1-alpha)) 
  
  #now create Z R&D resource use vector
  Z = matrix(0,T+1,1)
  Z[1:T] = gA[2:(T+1)]^gamma
  Z[1:T] = Z[1:T] * A[1:T]^(1+ ((1-delta)/rho) )
  Z[1:T] = Z[1:T] / ((gamma-1)*eta +1)

  #now create Northern and Southern consumption vectors
  C = matrix(0,T+1,1); Cstar = C
  C[1:T] = Y[1:T] - Z[1:T] - M[2:(T+1)] * (xM[2:(T+1)] + xMstar[2:(T+1)]) - R[2:(T+1)]*xR[2:(T+1)]
  Cstar[1:T] = Ystar[1:T] - I[2:(T+1)]*(xI[2:(T+1)]+xIstar[2:(T+1)]) - R[2:(T+1)]*xRstar[2:(T+1)]
  c=C/H; cstar = Cstar/Hstar;
  
  #now form the Northern HH Euler equation
  LHSNEULER = matrix(0,T-1,1); RHSNEULER = LHSNEULER
  LHSNEULER = (1+r[2:T])
  RHSNEULER = (1/beta) * (1+gH) * (c[2:T]/c[1:(T-1)])^sigma

  #now form the Southern HH Euler equation
  LHSSEULER = matrix(0,T-1,1); RHSSEULER = LHSSEULER
  LHSSEULER = (1+rstar[2:T])
  RHSSEULER = (1/beta) * (1+gH) * (cstar[2:T]/cstar[1:(T-1)])^sigma
  
  #now form the BT condition
  LHSBT = matrix(0,T-1,1); RHSBT = LHSBT
  LHSBT = q[2:T]
  RHSBT = ((phivec[2:T]/gA[2:T])*(1/human.cap.ratio))^(alpha/(2-alpha))
  RHSBT = RHSBT * ((1+r[2:T])/(1+rstar[2:T]))^((1-alpha)/(2-alpha))
  RHSBT = RHSBT * PSI

  obj = ((LHSNEULER - RHSNEULER)/RHSNEULER)^2
  obj = obj + ((LHSSEULER - RHSSEULER)/RHSSEULER)^2
  obj = obj + ((LHSBT - RHSBT)/RHSBT)^2
  obj=sum(obj)
   
  cat('error = ',obj,'\n')
  if (evaldum==0){return(obj)}
  if (evaldum==1){
    returnlist = list(gA,q,r,rstar,Y,Ystar,C,Cstar,c,cstar,H,Hstar,A)
    return(returnlist)
  }
  
}


c.fn.transition=cmpfun(fn.transition)

#now call the fn.transition function
pricevec = matrix(0,3*(T-1),1)
pricevec[1:(T-1)]=(q.phi+q.phi.prime)/2
pricevec[T:(3*(T-1))]=rBGP

paramvec = c(sigma,alpha,rho,eta,
             human.cap.ratio,beta,
             A1,H1,Hstar1,gH,
             gA,delta,T,rBGP,q.phi,q.phi.prime,
             phi,phi.prime)

lb = matrix(0,3*(T-1),1)
ub = matrix(0,3*(T-1),1)

#lb and ub for q
lb = pricevec*0.01
ub = pricevec*2.75

cat('Now doing the transition path optimization.')
transition=psoptim(pricevec,c.fn.transition,lower=lb,upper=ub,control=list(abstol=errortol,maxit=10000000),paramvec=paramvec,evaldum=0)

#debug(c.fn.transition)
returnlist = c.fn.transition(transition$par,paramvec,evaldum=1)

gAvec = returnlist[[1]]
qvec = returnlist[[2]]
rvec = returnlist[[3]]
rstarvec = returnlist[[4]]
Yvec = returnlist[[5]]
Ystarvec = returnlist[[6]]
Cvec = returnlist[[7]]
Cstarvec = returnlist[[8]]
cvec = returnlist[[9]]
cstarvec = returnlist[[10]]
Hvec = returnlist[[11]]
Hstarvec = returnlist[[12]]
Avec = returnlist[[13]]

plot(gAvec)
plot(Yvec[2:(T+1)]/Yvec[1:T]-1)
plot(Ystarvec[2:(T+1)]/Ystarvec[1:T]-1)
plot(cvec[2:(T)]/cvec[1:(T-1)]-1)
plot(cstarvec[2:(T)]/cstarvec[1:(T-1)]-1)
plot(1+rvec)
plot(1+rstarvec)
plot(Hvec/(Avec^((1-delta)/rho)))
plot(qvec)
plot(qvec*(1+rstarvec)-(1+rvec))

#now do no shock version
pricevec.noshock = matrix(0,3*(T-1),1)
pricevec.noshock[1:(T-1)]=q.phi
pricevec.noshock[T:(3*(T-1))]=rBGP

paramvec.noshock = c(sigma,alpha,rho,eta,
             human.cap.ratio,beta,
             A1,H1,Hstar1,gH,
             gA,delta,T,rBGP,q.phi,q.phi,
             phi,phi)

transition.noshock = psoptim(pricevec.noshock,c.fn.transition,lower=lb,upper=ub,control=list(abstol=errortol,maxit=10000000),paramvec=paramvec.noshock,evaldum=0)
returnlist.noshock = c.fn.transition(transition.noshock$par,paramvec.noshock,evaldum=1)

gAvec.noshock = returnlist.noshock[[1]]
qvec.noshock = returnlist.noshock[[2]]
rvec.noshock = returnlist.noshock[[3]]
rstarvec.noshock = returnlist.noshock[[4]]
Yvec.noshock = returnlist.noshock[[5]]
Ystarvec.noshock = returnlist.noshock[[6]]
Cvec.noshock = returnlist.noshock[[7]]
Cstarvec.noshock = returnlist.noshock[[8]]
cvec.noshock = returnlist.noshock[[9]]
cstarvec.noshock = returnlist.noshock[[10]]
Hvec.noshock = returnlist.noshock[[11]]
Hstarvec.noshock = returnlist.noshock[[12]]
Avec.noshock = returnlist.noshock[[13]]

#now, do the welfare comparisons
multiples = 10
Ttot = T*multiples
Wmat = matrix(0,Ttot,4)
Wmat[1:T,1]=cvec[1:T]
Wmat[1:T,2]=cstarvec[1:T]
Wmat[1:T,3]=cvec.noshock[1:T]
Wmat[1:T,4]=cstarvec.noshock[1:T]

betas=matrix(0,Ttot,1)
betas[1,]=1

for (t in (2:Ttot)){
  
  betas[t,]=betas[t-1,]*beta
  if (t>T){
    Wmat[t,] = Wmat[t-1,] * (1+gA)
  }
  
}

Wvec=c(sum(betas*log(Wmat[,1])),sum(betas*log(Wmat[,2])),sum(betas*log(Wmat[,3])),sum(betas*log(Wmat[,4])))

x = 100*(exp((1-beta)*(Wvec[1]-Wvec[3]))-1)
xstar = 100*(exp((1-beta)*(Wvec[2]-Wvec[4]))-1)


if (doplots==1){
  
  numper=25
  lwidnum=2  
  growth.lim=c(1.9,2.9)
  xrate.lim=c(0.4,0.75)
  r.lim = c(5,5.65)
  
  Ypercap=Yvec/Hvec
  Ygrowth = matrix(0,numper,1)
  Ygrowth[1] = 100*((1+gA)^(1/numyrs)-1)
  Ygrowth[2:numper]=100*((Ypercap[2:numper]/Ypercap[1:(numper-1)])^(1/numyrs)-1)
  
  Ystarpercap=Ystarvec/Hstarvec
  Ystargrowth = matrix(0,numper,1)
  Ystargrowth[1] = 100*((1+gA)^(1/numyrs)-1)
  Ystargrowth[2:numper]=100*((Ystarpercap[2:numper]/Ystarpercap[1:(numper-1)])^(1/numyrs)-1)
  

  pdf(file="semi_trans4.pdf")
  par(mfrow=c(2,2))
  
  plot(1:numper,100*((1+gAvec[1:numper])^(1/numyrs)-1),type="l",main="Variety Growth",xlab="Period",ylab="Annual %",col=1,lwd=lwidnum,ylim=growth.lim,xaxt="n")
  points(1:numper,100*((1+gAvec[1:numper])^(1/numyrs)-1),pch=1)
  lines(1:numper,100*((1+gA)^(1/numyrs)-1)*ones(numper,1),col=4,lwd=lwidnum)
  axis(1,at=c(1,5,10,15,20,25))
  
  plot(1:numper,qvec[1:numper],type="l",main="Southern Terms of Trade",xlab="Period",ylab="North/South",col=1,lwd=lwidnum,ylim=xrate.lim,xaxt="n")
  points(1:numper,qvec[1:numper],pch=1)
  lines(1:numper,q.phi.prime*ones(numper,1),col=4,lwd=lwidnum)
  lines(1:numper,q.phi*ones(numper,1),col=2,lty=2,lwd=lwidnum)
  axis(1,at=c(1,5,10,15,20,25))
  
  plot(1:numper,Ygrowth,type="l",main="Northern Output Growth",xlab="Period",ylab="Per Capita, Annual %",col=1,lwd=lwidnum,ylim=growth.lim,xaxt="n")
  points(1:numper,Ygrowth,pch=1)
  lines(1:numper,100*((1+gA)^(1/numyrs)-1)*ones(numper,1),col=4,lwd=lwidnum)
  axis(1,at=c(1,5,10,15,20,25))
  
  plot(1:numper,Ystargrowth,type="l",main="Southern Output Growth",xlab="Period",ylab="Per Capita, Annual %",col=1,lwd=lwidnum,ylim=growth.lim,xaxt="n")
  points(1:numper,Ystargrowth,pch=1)
  lines(1:numper,100*((1+gA)^(1/numyrs)-1)*ones(numper,1),col=4,lwd=lwidnum)
  axis(1,at=c(1,5,10,15,20,25))
  
  dev.off()
  
  numper=T
  Ypercap=Yvec/Hvec
  Ygrowth = matrix(0,numper,1)
  Ygrowth[1] = 100*((1+gA)^(1/numyrs)-1)
  Ygrowth[2:numper]=100*((Ypercap[2:numper]/Ypercap[1:(numper-1)])^(1/numyrs)-1)
  
  Ystarpercap=Ystarvec/Hstarvec
  Ystargrowth = matrix(0,numper,1)
  Ystargrowth[1] = 100*((1+gA)^(1/numyrs)-1)
  Ystargrowth[2:numper]=100*((Ystarpercap[2:numper]/Ystarpercap[1:(numper-1)])^(1/numyrs)-1)
  
  pdf(file="semi_trans6.pdf")
  par(mfrow=c(3,2))
  
  plot(1:numper,100*((1+gAvec[1:numper])^(1/numyrs)-1),type="l",main="Variety Growth",xlab="Period",ylab="Annual %",col=1,lwd=lwidnum,ylim=growth.lim)
  points(1:numper,100*((1+gAvec[1:numper])^(1/numyrs)-1),pch=1)
  lines(1:numper,100*((1+gA)^(1/numyrs)-1)*ones(numper,1),col=4,lwd=lwidnum)
  
  plot(1:numper,qvec[1:numper],type="l",main="Southern Terms of Trade",xlab="Period",ylab="North/South",col=1,lwd=lwidnum,ylim=xrate.lim)
  points(1:numper,qvec[1:numper],pch=1)
  lines(1:numper,q.phi.prime*ones(numper,1),col=4,lwd=lwidnum)
  lines(1:numper,q.phi*ones(numper,1),col=2,lwd=lwidnum)
  
  plot(1:numper,Ygrowth,type="l",main="Northern Output Growth",xlab="Period",ylab="Per Capita, Annual %",col=1,lwd=lwidnum,ylim=growth.lim)
  points(1:numper,Ygrowth,pch=1)
  lines(1:numper,100*((1+gA)^(1/numyrs)-1)*ones(numper,1),col=4,lwd=lwidnum)
    
  plot(1:numper,Ystargrowth,type="l",main="Southern Output Growth",xlab="Period",ylab="Per Capita, Annual %",col=1,lwd=lwidnum,ylim=growth.lim)
  points(1:numper,Ystargrowth,pch=1)
  lines(1:numper,100*((1+gA)^(1/numyrs)-1)*ones(numper,1),col=4,lwd=lwidnum)
  
  plot(1:numper,100*((1+rvec[1:numper])^(1/numyrs)-1),type="l",main="Northern Interest Rate",xlab="Period",ylab="Annual %",col=1,lwd=lwidnum,ylim=r.lim)
  points(1:numper,100*((1+rvec[1:numper])^(1/numyrs)-1),pch=1)
  lines(1:numper,100*((1+rBGP)^(1/numyrs)-1)*ones(numper,1),col=4,lwd=lwidnum)
  
  
  plot(1:numper,100*((1+rstarvec[1:numper])^(1/numyrs)-1),type="l",main="Southern Interest Rate",xlab="Period",ylab="Annual %",col=1,lwd=lwidnum,ylim=r.lim)
  points(1:numper,100*((1+rstarvec[1:numper])^(1/numyrs)-1),pch=1)
  lines(1:numper,100*((1+rBGP)^(1/numyrs)-1)*ones(numper,1),col=4,lwd=lwidnum)
  
  dev.off()
  
  
  
}


maxgA = gAvec[2]
halflife = which(gAvec<(gA + (maxgA-gA)/2))[2]


sink(file="summary.txt")
cat('Summary for Semi-Endogenous Version','\n')
cat('Balanced Growth Path Variety Growth, Annual % = ',100*((1+gA)^(1/numyrs)-1),'\n')
cat('Balanced Growth Path Interest Rate, Annual % = ',100*((1+rBGP)^(1/numyrs)-1),'\n')
cat('Balanced Growth Path I/Y(phi), % = ',100*IoverY.preshock,'\n')
cat('Balanced Growth Path I/Y(phi.prime), % = ',100*IoverY.postshock,'\n')
cat('phi, % = ',100*phi,'\n')
cat('phi.prime, % = ',100*phi.prime,'\n')
cat('Balanced Growth Path q(phi) = ',q.phi,'\n')
cat('Balanced Growth Path q(phi) = ',q.phi.prime,'\n')
cat('Northern Welfare Gains, % = ',x,'\n')
cat('Southern Welfare Gains, % = ',xstar,'\n')
cat('Halflife is ',halflife,' periods.','\n')
cat('Max gA ',100*((1+maxgA)^(1/numyrs)-1),'\n')
cat('Max impact (Max gA)-(BGP gA) = ',100*((1+maxgA)^(1/numyrs)-1)-100*((1+gA)^(1/numyrs)-1),'\n')

sink(file=NULL)
sink(file=NULL)

#save all of the results
save(file="semiendog.rda",list=ls())
end.time = Sys.time()
tot.time = end.time-start.time
cat('Total time was ',tot.time,' seconds.','\n')