%TRANSITION  Computes the transition path for
%
%   max  sum    beta^t (log c[i,t] + psi log (1-l[i,t])) N[i,t]
%      t  =0:inf
%
%   subject to
%
%     C[i,t] + Xk[i,t] + Xm[i,t] + NX[i,t]  <   = Y[i,t]
%     Y[i,t]      = A[i,t] [N[i,t] hatM[i,t]]^phi
%                       *[K[i,t]^alpha L[i,t]^(1-alpha)]^{1-phi}
%     K[i,t+1]    = (1-deltk) K[i,t] + Xk[i,t]
%     M[i,t+1]    = (1-deltm) M[i,t] + Xm[i,t]
%     hatM[i,t]   = M[i,t] + ome[i,t] M[-i,t]
%     NX[i,t]     = sum_{j\=i} (r[i,j,t] M[j,t] - r[j,i,t] M[i,t])
%     M[i,t]   >  = 0  (zeta is penalty parameter)
%     N[i,t]      = N[i,0]*(1+gamN)^t
%     A[i,t]      = (1+gamA)^t,  gamA=func(gamY,gamN)
%
%   and in equilibrium,
%     r[j,i,t]    = phi*ome[j,t] Y[j,t]/hatM[j,t]
%
%   Inputs:  I       number of countries 
%            alpha   capital share for non-technology input
%            beta    discount factor
%            deltk   depreciation rate of non-technology capital
%            deltm   dpereciation rate of technology capital
%            gamN    growth rate of N[i,t], all i
%            gamY    growth rate of Y[i,t], all i
%            phi     share of technology capital in production
%            psi     utility parameter
%            N       Ix1 vector of initial populations
%            T       number of periods
%            A       T*I x 1 vector of TFPs
%            sigma   T*I x 1 vector of openness parameters
%

%   Ellen McGrattan, 3-13-07

%
% Examples shown in the paper with alternative choices for N, A, sigma
% (Users can add their own examples)
%
disp('Transitions of Two Countries')
disp('  1.  Big is initially open and small joins')
disp('  2.  Two big countries, open at different times')
disp('  3.  Big is initially open and small joins and gains in A')
disp('  4.  Big is initially open and small joins without A gains')
disp(' ')
Case     = input(' Which case? ');
T        = 100;
if Case==1;
  %
  % Big country is open and small country joins
  %
  NI       = 10;
  NJ       = 1;
  ftanh    = [.995-tanh(-3.5:.2:-3.5+.2*(T-1))'*.005];
  xtanh    = (1-ftanh)*99;
  sigmat   = [xtanh(T)*ones(T,1),xtanh];
  At       = ones(T,2);
end;
if Case==2;
  %
  % Two big countries, country 2 opens second
  %
  NI       = 10;
  NJ       = 10;
  ftanh    = [.995-tanh(-1:.08:-1+.08*(T-1))'*.005];
  ftanh2   = [.995-tanh(-1.5:.07:-1.5+.07*(T-1))'*.005];
  sigmat   = [.45*(1+(1-ftanh)*105),.17*(1+(1-ftanh2)*431)];
  At       = ones(T,2);
end;
if Case==3;
  %
  % Big country is open and small country joins and has A gains
  %
  NI       = 10;
  NJ       = 1;
  ftanh    = [.995-tanh(-3.5:.2:-3.5+.2*(T-1))'*.005];
  xtanh    = (1-ftanh)*99;
  sigmat   = [xtanh(T)*ones(T,1),xtanh];
  At       = [ones(T,1),.9+.1*xtanh];
end;
if Case==4;
  %
  % Big country is open and small country joins and has no A gains
  %
  NI       = 10;
  NJ       = 1;
  ftanh    = [.995-tanh(-3.5:.2:-3.5+.2*(T-1))'*.005];
  xtanh    = (1-ftanh)*99;
  sigmat   = [xtanh(T)*ones(T,1),xtanh];
  At       = [ones(T,1),.9*ones(T,1)];
end;
%
% Fixed parameters
%

alpha    = .3;
beta     = .98;
deltk    = .053;
deltm    = .1;
gamN     = .01;
gamY     = .029392;
phi      = .07;
psi      = 2.5;
constant = [alpha;beta;deltk;deltm;gamN;gamY;phi;psi];
AI       = At(1,1);
AJ       = At(1,2);
VI       = balgrowth(NI,AI,sigmat(1,1),constant);
VJ       = balgrowth(NJ,AJ,sigmat(1,2),constant);

%
% Starting point for transition and guess of solution
%
k0       = [VI(6);VJ(6)];
m0       = [VI(8);VJ(8)];
y0       = [VI(1);VJ(1)];
c0       = [VI(2);VJ(2)];
xk0      = [VI(3);VJ(3)];
xm0      = [VI(4);VJ(4)];
l0       = [VI(5);VJ(5)];

for i=1:T;
  Vtem   = balgrowth([NI;NJ],At(i,:)',sigmat(i,:)',constant);  
  c(i,:) = Vtem(:,2)';
  l(i,:) = Vtem(:,5)';
  xk(i,:)= Vtem(:,3)';
  xm(i,:)= Vtem(:,4)';
end;

%Vtem    = balgrowth([NI;NJ],At(1,:)',sigmat(1,:)',constant);  % Alternate
%k0      = Vtem(:,6);                                          % initial
%m0      = Vtem(:,8);                                          % condition

ikt      = find(xk>0);
imt      = find(xm>0);
param    = [2;alpha;beta;deltk;deltm;gamN;gamY;phi;psi;0;NI;NJ;k0;m0;T;
            At(:);sigmat(:);length(ikt);length(imt);ikt;imt];
Xt       = [c,l];
Xt       = [Xt(:);xk(ikt);xm(imt)];
Xt       = secant('transfoc',Xt,param);
[r,Vt]   = transfoc(Xt,param);
z        = zeros(T,2);
yt       = z; yt(:)  = Vt(:,1);
ct       = z; ct(:)  = Vt(:,2);
xkt      = z; xkt(:) = Vt(:,3);
xmt      = z; xmt(:) = Vt(:,4);
lt       = z; lt(:)  = Vt(:,5);
kt       = z; kt(:)  = Vt(:,6);
mhat     = z; mhat(:)= Vt(:,7);
mt       = z; mt(:)  = Vt(:,8);
nxt      = z; nxt(:) = Vt(:,9);
knt      = [kt(2:T,:);kt(T,:)];
mnt      = [mt(2:T,:);mt(T,:)];

if any(xkt(:)<0 | xmt(:)<0);
  newnegk  = find(xkt<0);
  newnegm  = find(xmt<0);
  Xtchk    = Xt;
  newikt   = find(xkt>0);
  newimt   = find(xmt>0);
  parchk   = [2;alpha;beta;deltk;deltm;gamN;gamY;phi;psi;0;NI;NJ;k0;m0;T;
              At(:);sigmat(:);length(newikt);length(newimt);newikt;newimt];
  Xtchk    = [Xt(1:4*T);xkt(newikt);xmt(newimt)];
  Xtchk    = secant('transfoc',Xtchk,parchk);
  [r,Vt]   = transfoc(Xtchk,parchk);
  z        = zeros(T,2);
  yt       = z; yt(:)  = Vt(:,1);
  ct       = z; ct(:)  = Vt(:,2);
  xkt      = z; xkt(:) = Vt(:,3);
  xmt      = z; xmt(:) = Vt(:,4);
  lt       = z; lt(:)  = Vt(:,5);
  kt       = z; kt(:)  = Vt(:,6);
  mhat     = z; mhat(:)= Vt(:,7);
  mt       = z; mt(:)  = Vt(:,8);
  nxt      = z; nxt(:) = Vt(:,9);
  knt      = [kt(2:T,:);kt(T,:)];
  mnt      = [mt(2:T,:);mt(T,:)];
end;
                                                                                
t        = 1:T;
figure(1)
set(gca,'Fontsize',14)
plot(t,sigmat(:,1),'Linewidth',2,'Color','b')
line(t,sigmat(:,2),'Linewidth',2,'Color','r')
title('Inputs: Degree of Openness \sigma_i_t (in TFP space)','Fontsize',14)
figure(2)
set(gca,'Fontsize',14)
plot(t,sigmat(:,1).^(1/phi),'Linewidth',2,'Color','b')
line(t,sigmat(:,2).^(1/phi),'Linewidth',2,'Color','r')
title('Inputs: Degree of Openness \omega_i_t (in fraction-let-in space)', ...
      'Fontsize',14)
figure(3)
set(gca,'Fontsize',14)
plot(t,At(:,1),'Linewidth',2,'Color','b')
line(t,At(:,2),'Linewidth',2,'Color','r')
title('Inputs: Total Factor Productivity A_i_t','Fontsize',14)
figure(4)
set(gca,'Fontsize',14)
plot(t,ct(:,1)/c0(2)*100,'Linewidth',2,'Color','b')
line(t,ct(:,2)/c0(2)*100,'Linewidth',2,'Color','r')
title('Per Capita Consumption Relative to c_2_,_0(1+\gamma_y)^t', ...
      'Fontsize',14)
ylabel('%')
figure(5)
set(gca,'Fontsize',14)
plot(t,yt(:,1)./lt(:,1)*l0(2)/y0(2)*100,'Linewidth',2,'Color','b')
line(t,yt(:,2)./lt(:,2)*l0(2)/y0(2)*100,'Linewidth',2,'Color','r')
title('Productivity Relative to y_2_,_0/l_2_,_0(1+\gamma_y)^t','Fontsize',14)
ylabel('%')
figure(6)
set(gca,'Fontsize',14)
ym0 = (y0-xm0./[NI;NJ])./l0;
YMt = [(yt(:,1)-xmt(:,1)/NI)./lt(:,1),(yt(:,2)-xmt(:,2)/NJ)./lt(:,2)];
plot(t,YMt(:,1)/ym0(2)*100,'Linewidth',2,'Color','b')
line(t,YMt(:,2)/ym0(2)*100,'Linewidth',2,'Color','r')
title(['Measured Productivity Relative to (Y_2_,_0-X_m_2_,_0)/L_2_,_0', ...
       '(1+\gamma_Y)^t'],'Fontsize',14)
ylabel('%')
figure(7)
set(gca,'Fontsize',14)
plot(t,mt(:,1)/(NJ*y0(2)),'Linewidth',2,'Color','b')
line(t,mt(:,2)/(NJ*y0(2)),'Linewidth',2,'Color','r')
title('Technology Capital Relative to Y_2_,_0(1+\gamma_Y)^t','Fontsize',14)

