% BIGvSmall   computes the steady states for the 2-country case
%             and generates Figure 8 and 9 of 
%
%                  Openness, Technology Capital, and Development
%                  by Ellen McGrattan and Edward Prescott
%

%             Ellen McGrattan, 3-15-07
%             Revised, ERM, 9-22-07

%
% Computing the steady state for 2-country case
%

A        = input(' Input the 2x1 A vector: ');
N        = input(' Input the 2x1 N vector: ');
A1       = A(1);
A2       = A(2);
N1       = N(1);
N2       = N(2);
alpha    = .3;
beta     = .98;
deltk    = .053;
deltm    = .1;
gamN     = .01;
gamY     = .029392;
phi      = .07;
Z1       = A1^(1/(1-alpha*(1-phi)));
Z2       = A2^(1/(1-alpha*(1-phi)));
super1   = (1-alpha*(1-phi))/((1-alpha)*(1-phi));
super2   = phi/(1-alpha*(1-phi));
super3   = super1*super2;
gamy     = (1+gamY)/(1+gamN)-1;
rho      = (1+gamy)/beta-1;
ak       = alpha*(1-phi)/(rho+deltk);
am       = phi/(rho+deltm);
ycoef    = (ak^(alpha*(1-phi))*am^phi)^(1/((1-alpha)*(1-phi)));

%
% Critical value of omega
%
omecrit  = (Z2*N2/(Z1*N1))^super1;
ome      = linspace(0,1,500)';
i1       = find(ome<=omecrit);
i2       = find(ome>omecrit);
ome1     = ome(i1);
ome2     = ome(i2);
tem      = (Z1*N1+Z2*N2*ome2.^super2).^super3;
YdivZN   = ycoef*[[(1+ome1)*[Z1*N1,Z2*N2]].^super3;
                  [tem,tem.*ome2.^super2]];
Y1       = YdivZN(:,1)*Z1*N1;
Y2       = YdivZN(:,2)*Z2*N2;
M1       = phi/(rho+deltm)*[(Y1(i1)-ome1.*Y2(i1))./(1-ome1);
                           [Y1(i2)+Y2(i2)]];
M2       = phi/(rho+deltm)*[(Y2(i1)-ome1.*Y1(i1))./(1-ome1);
                           0*ome2];
C1       = [Y1(i1)*(1-deltk*ak)-deltm*am*(Y1(i1)-ome1.*Y2(i1))./(1-ome1)- ...
            phi*ome1.*(Y1(i1).*M2(i1)./(M1(i1)+ome1.*M2(i1))- ...
            Y2(i1).*M1(i1)./(M2(i1)+ome1.*M1(i1)));
           (1-deltm*am)*(Y1(i2)+Y2(i2))-(1-phi)*Y2(i2)-deltk*ak*Y1(i2)];
C2       = [Y2(i1)*(1-deltk*ak)-deltm*am*(Y2(i1)-ome1.*Y1(i1))./(1-ome1)- ...
            phi*ome1.*(Y2(i1).*M1(i1)./(M2(i1)+ome1.*M1(i1))- ...
            Y1(i1).*M2(i1)./(M1(i1)+ome1.*M2(i1)));
            (1-phi-deltk*ak)*Y2(i2)];
Xk1      = [deltk*ak*Y1(i1);
            deltk*ak*Y1(i2)];
Xk2      = [deltk*ak*Y2(i1);
            deltk*ak*Y2(i2)];
Xm1      = deltm*M1;
Xm2      = deltm*M2;
NX1      = [phi*ome1.*(Y1(i1).*M2(i1)./(M1(i1)+ome1.*M2(i1))- ...
            Y2(i1).*M1(i1)./(M2(i1)+ome1.*M1(i1)));
            -phi*Y2(i2)];
NX2      = [phi*ome1.*(Y2(i1).*M1(i1)./(M2(i1)+ome1.*M1(i1))- ...
            Y1(i1).*M2(i1)./(M1(i1)+ome1.*M2(i1)));
            phi*Y2(i2)];
ome      = [ome1;ome2];
figure(1)
set(gca,'Fontsize',14)
plot(ome,YdivZN(:,1)/YdivZN(1,2)*100,'-', ...
     ome,YdivZN(:,2)/YdivZN(1,2)*100,'-')
title('Indices of Output per Effective Person Relative to Trend', ...
      'Fontsize',14)

figure(2)
set(gca,'Fontsize',14)
plot(ome,M1./Y1,'-',ome,M2./Y2,'-')
title('Technology Capital Relative to Output', ...
      'Fontsize',14)

fig8 = [ome, YdivZN(:,1)/YdivZN(1,2)*100, YdivZN(:,2)/YdivZN(1,2)*100];
fig9 = [ome, M1./Y1, M2./Y2];
